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SUMMARY 


The instability of the fabrication process of optical fibers introduces 
both ellipticies and stress anisotropy. These perturbations are the causes 
for birefringence in single-aode optical fibers which have been researched 
extensively. In this research, the node propagations in optical fibers with 
anisotropic elliptical core have been Investigated. 

The exact characteristic equation for anisotropic elliptical optical 
fibers can be obtained for odd and even hybrid Bodes in terns of infinite 
deterainants utilizing Mathleu and aodlfied Mathleu functions. The exact 
characteristic equation is applicable to elliptical fibers that have any 
elllpticity. A siapllfied characteristic equation can then be obtained by 
applying the veakly guiding approxination such that the difference in the 
refractive indices of the core and the cladding is snail. It has been shovn 
that significant sinplif icatlon can be achieved under this approximation. 

The siapllfied characteristic equation is used to compute the 
propagation constants for the anisotropic elliptical fiber. The expression 
for the pover carried by the fiber is also obtained and numerical results are 
presented. These results may be used to approximate a number of different 
shapes of fibers. 


vi 



1. INTRODUCTION 


The circular optical fiber is one of the Boat studied media for long 
distance coaaunicatlon. An optical fiber consists of a core of a dielectric 
material in vhich the refractive index is higher than the refractive Index of 
the cladding. However, a cladded fiber is often tlaes aodeled as a dielectric 
rod when the cladding radius is large enough such that the guided aode fields 
will decay to insignificant values at the outer boundary of cladding. The 
theory of optical fibers of this type is veil understood and has been 
described in detail in the previously published research! 1,2] . 

In an effort to obtain a lov-loss fiber, Honerle!3] carried out an 
experiaental study of doubly clad fibers in vhich the refractive index of 
inner cladding vas less than that of core and outer cladding. This study 
shovs that the optlaua doping levels in the core of doubly clad fibers are 
less than those required by dispersion-free singly clad fibers. This leads to 
a smaller propagation loss, since the scattering losses decreased vith a 
decreasing dopant concentration in fibers. 

It is interesting to note, hovever, that the Instability in the 
fabrication process nay introduce ellipticities in the optical fibers. This 
lack of circular synnetry is one of the causes for birefringence in single- 
node optical fibers; such a birefringent fiber is also called a single- 
polarization single-node fiber! 4]. These birefringent fibers are inportant 
for systems utilizing such fibers as fiber optic sensors and for predicting 
the transmission bandvldth reduction caused by group-delay differences between 
orthogonally polarized nodes. 

The birefringence due to elllptlclty has been studied 
experimentally! 5, 6] and the neasured data have been compared with those 
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equation for an unlaxlally anisotropic circular rod for hybrid modes of 
excitation. Analytic solutions for the circular fiber vhen both core and 
cladding consist of uniaxial aaterlal was presented by Tonnlng(32]. This 
study Indicates that the cut-off frequency for the lovest-order mode is not 
affected by the cladding anisotropy. 

Vhen the circular cross-section of the fiber is deformed into a 
nonclrcularly syaaetric profile, a single aode In a circular fiber may split 
into two aodes vith different polarizations and propagation velocities! 33] . 
This has been experimentally verified by employing the near-field aethod(3U 
and the spectral polarization method!35). Cozen and Dyott!36] obtained the 
cut-off frequency of the first higher order aode in an elliptical fiber fron 
an approximate characteristic equation. However, the limitation of their 
results is described by Clterne(37] and Rengarajan{38) . The cut-off 
characteristic has also been obtained by solving the exact characteristic 
equation in terms of Mathleu functions and modified Mathleu functions (38, 39 ] 
and by applying the mode-matching method(27J or the critical wavelength shift 
foraular method (401. However, there exists a disagreement in the 
interpretation of their results, especially in the region where the ratio 
between minor axis and major axis is small; Saad(41] presented possible 
reasons for these differences. 

Por most of the practical fibers used as optical communication lines, 
the simplification of the characteristic equation is possible by applying the 
weakly guiding approximation such that the difference in the refractive 
Indices of the core and the cladding is small!26,33,42,43] . It is shown that 
the error introduced by this simplification is less than 10% even vhen the 
difference in the refractive indices is equal to two. The perturbation method 
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can also be applied to study the polarization effects in aulti-aoded fibers 
when the fiber is weakly guiding and/or weakly anisotropic t 44] . It is also 
possible for aulti-aode fibers to siaplify the characteristic equation by 
applying a perturbation aethod based on the far-froa-cut-of f approxiaatlon as 
shown by Paul and Shevgaonkar(45] in a study of circular fiber with uniaxial 
anisotropy. This approxiaatlon is useful for the aulti-aode propagation in 
optical waveguides since the lover order nodes that carry aost of the power 
could be considered in the far-from-cut-off region. 

For the aore general case of biaxial anisotropic waveguide, an analytic 
solution of the field equations is not possible even for waveguides with 
sinple geometries . However, the fields and propagation constants can be 
obtained by applying the numerical techniques discussed above. The 
propagation constants can also be coaputed by using a coupled node theory. 

This coupled node approch has been applied for the study of node propagation 
in rectangular guides(46] and cylindrical fibers(47]. 

is it has been shown through the previous discussion, the instability of 
the fabrication process of optical fibers Introduces both elllptlcities and 
stress anisotropy. Also, the results obtained for an elliptical optical fiber 
aay be used to approximate a number of different shapes of fibers. It can 
take the shape of a circular fiber or that of a flat tape fiber depending upon 
the eccentricity of the elliptical fiber. Hence, it is proposed to 
investigate the mode propagation in elliptical optical fibers containing 
uniaxial anisotropic media. In this study, the fiber will be modeled as a 
dielectric elliptical rod, since the departure of the cladding's cross-section 
froa circular form can be ignored In the case of large dlaension of cladding 
radius. The exact characteristic equation for the anisotropic elliptical 
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fiber having any elllptlclty vill be obtained using the series of Mathleu and 
■odlfied Mathleu functions. A simplified characteristic equation vill then be 
obtained by applying the veakly guiding approximation and the computed results 

vill be presented. 



2. WAVE EQUATION IN ELLEIPTICAL COORDINATES 
In solving Maxwell's equations, the wave equation In the waveguide can 
be expressed in the orthogonal curvilinear coordinates ( |?2» 48 

(l/fti 2 ) ( a 2 E z /aB! 2 ) ♦ U/£2 2 H » 2 E z /aJf2 2 > 

(2.1) + (l/ei«2)f <3( £ 2 / *1 ) /« SfL * E z /a»i ) + 

+ kl 2 E z « 0 

where kj is a constant and and 82 are Multiplying factors depending upon 
the particular coordinates [ 48 ] . is replaced by -10. An Identical 

equation can be obtained for H z . 

For the elliptical coordinates shown in Figure 1, 

(2.2) « f , 1 

and 

(2.3) ®1 * ^2 “ 9 t( cosh - cos 2? )/2J*. 

By substituting Eqs.(2.2) and (2.3) into Eq.(2.1), the following equation is 
obtained 

(2.4) 3 2 E z /ai 2 ♦ &E z /d f 2 ♦ 2k 2 ( cosh 2 £ - cos 2*? )E Z « 0 

with 2k * kjq. Then Eq.(2.4) is the two-dlaentional wave equation in 
elliptical coordinates. 

If we let the solution of Eq. (2.4) be E z ( £ , •? ) *'!'(%)♦<*?), Eq. (2.4) 
becoaes 

(2.5) fd 2 f/d| 2 *■ td 2 4/d* 2 + 2k 2 ( cosh 2 1 - cos 2V)'f'<t>* 0 
Dividing Eq.(2.5) by 4$ yields 

(2.6) (1/p d 2 Y/dfc 2 ♦ 2k 2 cosh 2l) + (1/4 d 2 4/d? 2 - 2k 2 cos 2fr) =0. 

Since the equations in the parenthis are independent to each other, we obtain 

(2.7) d 2 4/d , 7 2 + ( a - 2k 2 cos 2 f)$ s 0 
and 
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(2.8) dty/dt 2 - ( a - 2k 2 cosh 21)4*= 0 


E, 


where a la the separation constant. Then the solutions for Bq .(2.4) are 

(2.9) 

for k 2 > 0 

( 2 . 10 ) 


{ 


1 


Ce B ( £ , k 2 )ce B ( , k 2 ) 

(even) 

Se B ( | , k 2 )se B ( V , * 2 ) 

(odd) 

Fek B ( £ , k 2 )ce B ( *? , k 2 ) 

(even) 

Gek B ( | , k 2 )se B ( •? , k 2 ) 

(odd) 


) 


for k 2 < 0. 

Simllary, the solutions for H z can be obtained using the aethod discussed 


above . 



3. CHARACTERISTIC EQUATION 


The geometry shown in Figure 1 consists of an uniaxialy anisotropic 
elliptical rod with a permittivity tensor 


(3.1) 


E 1 


Cj 0 0 

0 0 

0 0 gC} 


which is embedded in a lossless dielectric medium of permittivity €o. 


anisotropic parameter g indicates the effect of anisotropic dielectric. 


The 

The 


anisotropic parameter is unity for isotropic case. 

It has been shown that in order to satisfy the boundary conditions 
completely both longitudinal electric and magnetic fields must be present, 
thus only hybrid type modes exist in elliptical fibers(2). Furthermore, due 
to the asymmetry of the elliptical cylinder, two types of modes exist and they 
are designated as an even type modes and an odd type modes. 


3.1 EVEN HOPES 

Assuming the t-z dependence of e^ w ^“® z ^for all field compornents, where 

0 is the propagation constant and w is the angular frequency, the axial 

components of the field for even modes are 

Eji * ? *lm Se m( ^ t Y le^ 8e mt ^ t Y lc 2 ^ 

(3.2) *"* 

Hji ^ t Y lh 2 ) ce mt *? , Y ih 2 ) 

for 0 i i s 6o 

E z2 ■ ?. t A 2m Ge)c m< Y 2 2 )Be B (?, Y 2 2 > 

H Z 2 ■ X k B2i|Pek m ( * ' Y 2 2jce m ( V# *2 2 > 

for lo S I < • 
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where *1. and B} B , 1 * 1,2 are arbitrary constants, and 
Y le 2 “ q 2 /4(v 2 )i g«i - gB 2 ) 

(3.4) t lh 2 ■ q 2 /4(v 2 M« 1 - 0 2 ) 

T 2 2 - q2/4(v 2 ye 0 - $2) 

q is the seaifocal length of the ellipse and M is the peraeabllity. 
The transverse field components are 
for 0 s i s to 

(3.5) Eh » ”i/ ( v 2 )i€ 1 - B 2 )L 

( & ( t 1 Y le 2 ) 8 *a^» Y le 2 ^ 

♦ X B l« Ce « ( 1 ' T lh 2 > c *a'( Y lh 2 H 

(3.6) Ei,! - -l/(vVi - 0 2 )L 

( B J^AimSeat E , Y ie 2 )*®a'^^* Y le 2 ^ 

" X B l« Ce »' ( £ ' Y lh 2 > c ®a<?» Y lh 2 >J 

(3.7) Hj! = -i/(v 2 H«! - 0 2 )L 

( *l«Se B ( E , Y le 2 )se B ' ( V , Y le 2 ) 

+ 8 ?' B la Ce M '( i, T lh 2 )ce B C?, T lh 2 )) 

(3.8) - -1/(v 2 M€i - B 2 )L 

I v«l5,il.Se B Ml, Yi e 2 )se B <7, T^ 2 ) 

+ ^ ' Y lh 2 ^ ce a’^» Y lh 2 ^ 

for € 0 s i < • 

(3.9) E| 2 ■ -i/(v 2 p€ 0 - B 2 )L 

( 0£A 2B Gek B Mi, T 2 2 )se B C? # t 2 2 ) 

♦ vn ? % B 2B Pek B ( t , r 2 2 )ce B *(7, y 2 2 )) 

(3.10) Eij 2 « -i/(v 2 n6o - B 2 )L 

( B i ? | i 2B Gek B ( i, Y 2 2 )se B * ( 7 , t 2 2 ) 

" V M Xe B 2B Pek a ,( Y 2 2 > ce a<?» y 2 2 H 
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(3.11) H| 2 * -1/(w 2 H€ 0 " 0 2 > l 

( -v« 0 j: A2a G « k m< E , Y 2 2 )se B '(7, Y 2 2 ) 

+ P? w B 2> Fek 1 , ( 6 , Y 2 2 )ce B (7# Y 2 2 )l 

(3.12) H$ 2 - -1/(w 2 M€ 0 - 0 2 )L 

( v€ 0| £ ( A 2B Gek B '( % , Y 2 2 )se B (7, y 2 2 ) 

♦ 0 E B 2a Fek B ( % , Y 2 2 )ce B '( 7 , Y 2 2 )J 

vhere 

L * ql( cosh2£ - cos27 )/2J* 

and the derivative vlth respect to £ or 7 is denoted by the prime. 

The boundary conditions require that the tangential B and H fields be 
continuous at the dielectric discontinuities. Equating the tangential fields 
at the boundary surface, f ■ £©/ 9 iveB 


(3.13) 

JS A lB Se B ( l 0 )se B ( 7 ) - X A 2 m G « k a< £ 0 >»«a*< V > 

(3.14) 

J? e B la c «m< io> ce mC?) B X B 2 » Fek s ( £o> c «a*< V > 

(3.15) 

l/(v 2 ^e 1 - 02 ) [ 0 j! ( *lm Se m{ lo^em't 1 ?* 

- vnX B U Ce » ,( g o> ce m< 12 >1 
- l/(v 2 pe 0 - 0 2 ) l 0 f S ( * 2 m Ge *i< £o)se B *'(7) 

- vnX B 2 m Fek m'( £o> ce m*< ? H 

(3.16) 

l/(v 2 H«l - 0 2 ) l vfi X a 1 b s *b' < fc 0 >sem<7) 

♦ 0 f ? o Bi B Ce > ( £ 0 )ce B ' (7)1 

- 1 / (v 2 p € 0 - 0 2 ) 1 vf 0 X* 2 a Gek m'< £o> 8 «a*<7 ) 

♦ B ^ # B 2 a Fek m^ £o^ ce a V ^ 

vhere the 

following abbreviations have been used. 

(3.17) 

Se B ( | 0 ^ * S*a( to# T le 2 ) 

(3.18) 

se B ( 7 ) * se B ( 7 , Y le 2 ) 

(3.19) 

Ce B ( f 0 ) * Ce B ( £q, Y^h 2 ) 

(3.20) 

ce B ( 7 ) * ce B ( 7 , ^lh 2 ) 
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(3.21) Gek a ( l 0 ) - Gek a ( l 0 , Y 2 2 ) 

(3.22) se B *( n ) - se a (»7, t 2 2 ) 

(3.23) Fek a ( % 0 ) « Fek a ( £ 0 , y 2 2 ) 

(3.24) ce B *( «? ) » ce B ( ^ , Y 2 2 ) . 

Multiplying both sides of Eqs. (3.13) and (3.16) by se n ( •? ) and Eqs. (3.14) 
and (3.15) by ce n ( *2 ), integrating with respect to V fro* 0 to 2x, and 
applying the orthogonality relations of the angular Mathleu functions 

(3.25) / ce B ce n • 0 if a 4 n 

e 

leads to 


(3.26) 

A ln Se n * *2m ^a^n 


(3.27) 

B ln Ce n * X B 2a p ®*m a m,n 


(3.28) 

B X B la Sa a yn /» - *HB ln Ce n ' - 



enh 2 /»2 2 X A 2« Gek » y n,« " 

"M Y ih 2 / Y 2 2 £. B 2a Pa *i ,a n,a 

(3.29) 

v€ l A ln Se n * B la Ce a^n,a * 



vt o Y lh 2 / Y 2 2 X 

♦ er lh 2 /Y 2 2 j£ B 2a p ®*a*n,a 


The priae over the summation sign is used to indicate that either odd or even 
values of a are used accordingly as to whether n is odd or even. 

«a,n* 0a,n/ ^a^n and y m,n ar * 9 iven by the following 

(3.30) « B/ n « £*ce B *($) ce n (S) d? / 4*ce n 2 (*) 

(3.31) 0 B/I) » J*se B # (*) se n (*) d* / d * 

(3.32) H» B/n «J?fce B '(*) se n <*) dl / 4*se n 2 m d* 

(3.33) y B/fl • £% e B • ( ■? ) ce n (9) d$ / d*. 

Making use of Eqs. (3.26) and (3.27), Eqs. (3.28) and (3.29) yields tvo 
sets of infinite homogeneous equations 

Jr,' A 2m s a,n + X B 2a t a / n " 0 
X A2aga » n + X B 2* h a,n ■ ° 


(3.34) 
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vhere 

(3-35) 9b, n * “Cl - nh 2 / Y 2 2 >°* k m< to> *m,r y r,n 

(3.36) h B ,n " I Ftkgl l^)Ce a '( |q)/CS|( £©) “ ^o^lh 2 ^^2 2 1 

(3.37) s B , n * v0 B/n /0 l e l 6ell »( £o> 8e m'< So>/ 8e m< M 

- € 0 Gek B ’( £o)W/ T 2 2 l 

(3.38) t 1#n - (1 " n h 2 A2 2 >Fek,( ^ 0> *■»* ^» n 

Por a nontrivial solution, the Infinite determinant of Eg. (3.34) must vanish. 
The propagation constant 0 can then be determined from the roots of this 
infinite determinant. 

The Infinite determinant for odd values of m and n is 

sll til s31 t31 - - 

gll hll g31 h31 - - 

sll tl3 s33 t33 - - 

gl3 hl3 g33 h33 - - 

: : s : - - 

: : : : 

and for even values of a and n is 

hOO g20 h20 g40 - - 

t02 s22 t22 s42 - - 

h02 g22 h22 g42 - - 

t04 s24 t24 s44 - - 

s j : s - - 
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3.2 ODD MODES 


The axial coaponents of the field for odd Bodes are 


(3.41) 


*zl T le 2)ce B (7, r le 2) 

«zl “J(, Bi B 3e B ( f , Y lh 2 ) se a( *1 • Y lh 2 ) 


for 0 s | S £q 


(3.42) 


*z2 “ 1 ?,*2» r * ,t B< t , Y 2 2 > ce a<^» Y 2 2 > 
«z2 -I B 2 «6ek.(£ # T 2 2 )ae.(? # Y 2 2 ) 


for l 0 s | < • 

vhere A} B and B^ n , 1 « 1,2 are arbitrary constants, and T^ e 2 , Yu, 2 , and y 2 2 
are given in Bq. (3.4) . 

The transverse field coaponents are 
for 0 i ? i Iq 

(3.43) Bn * -i/tv^Ci - 0 2 )L 

l fl£.*laC e.ME, Y le 2 )ce B (7, Y le 2 ) 

♦ w ^X B la Se » ( l ' Y lh 2 )®®a’<?. Y lh 2 > 1 

(3.44) B^i * “1/ ( v 2 p€^ - e 2 )L 

( 0 ?^*la c *a^ ^ • Y le 2 ) ce a' ( ^ / Y le 2 ^ 

“ ^,®la^ e a' ( f • Y lh 2 ^* e a^ • Y lh 2 ^ 

(3.45) Hfc • -i/wVi - B 2 )L 

1 ■ vC lX A i* Ce “ ( 1 ' Y le 2 ) c «a'< Y le 2 > 

♦ 0 Jf,B lB Se B '( £, Yi h 2 )se B e?, Wn 

(3.46) « -i/(v 2 H«i - 0 2 )L 

( J? e &i B Ce B ' ( E / Yi e 2 )ce B ( V , Y le 2 ) 

+ ®^i®la® e a^ • Y lh 2 )® e B'(^» Y lh 2 )l 



for £ 0 i I < • 

(3.47) E|2 ■ -l/(w2p€ 0 - 0 2 )L 

( P ^, 0 A 2m Fek «’< T2 2 )ce i |( ^ , *2 2 ) 

♦ Wi | ? ( B 2> Gek 1 ( f , Y 2 2 )»e B , (i, Y 2 2 )l 

(3.48) Eij2 ■ -l/lvfyCo - 0 2 )L 

l P j[ 0 *2M F « k s< * ' T 2 2 )ce B MV, t 2 2 ) 

- vpJ? ( B 2B Gek B '( t , V 2 2 )8e«( ^ $ *2 2 > ] 

(3.49) H|2 ■ -i/(v 2 Jl€ 0 - 0 2 )L 

l -v£ 0 j^A 2B Fek B ( £, T 2 2 )ce B '(?, T 2 2 ) 

♦ 0 < 5 i B2 B Gek 1 ' ( 1 1 Y2 2 )»e B (7# Y2 2 )l 

(3.50) H»^2 ■ -l/( v 2 H€ 0 “ P 2 > l 

( ^ 1 T 2 2 ) ce B ( *1 ' *2 2 ) 

♦ 0 | ? f B2 B Gek B ( t , T2 2 )ae B '(7f T 2 2)J 

The derivative with respect to £ or *? is denoted by the prise. 

Equating the tangential fields at the boundary surface, £ ■ £o» 9l yes 

(3.51) Jf e *l« c «a< &o) ce a< *? ) * Ji, A 2« Fek s ( &o) ce s ,(l ? 1 

(3.52) Jf,Bi B Se B ( i 0 >«e B (7) »X B 2» Gek * ( ^o> se » t( * 1 

(3.53) l/tv 2 )!^ - 0 2 ) I aJ|,A lB Ce B ( folce*'^) 

- vpJE B lB Se B '( £ 0 )se«( *1 >1 

* l/(w 2 )i€ 0 - 0 2 ) I 0 J? c A 2B Fek B ( £ 0 )ce B *'C7) 

- v)i | ? | B2 B Gek B , ( t 0 )se B *( •? )) 

(3.54) 1/(v 2 M € i - 0 2 ) ( v€i ? et Ai B Ce B ' ( £ 0 )ce B ( *1 ) 

♦ eJ?,Bi B Se B ( £ 0 ,se « ,( ‘? )1 

* l/(v 2 M«o * P 2 > 1 y €o^ A 2a Fek « ,( £o> ce » #( ^ } 

+ 0 S B 2B Gek B ( £ 0 )se B *' ( ^ ) 1 
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The 

! abbreviations 



(3. 

55) 

Ce B ( ) ■ 

Ce B ( Sor 

*le 2 ) 

(3. 

56) 

ce B ( *7 ) ■ 

ce B ( *7 , 

Tie 2 ) 

(3. 

57) 

8««(t 0 ) ■ 

S*i( t 0 t 

Tlh 2 ) 

(3. 

58) 

se B ( «? ) ■ 

»e B ( *7 1 

Tlh 2 ) 

(3. 

59) 

rek B (fo) 

« Pek B ( i 

0/ t 2 2 

(3. 

60) 

ce B *( •? ) 

• ce B ( *?# 

Y 2 2 ) 

(3. 

61) 

Oek B (l 0 ) 

- Oek B ( % 

0/ T 2 2 

(3. 

62) 

•e B # (2 ) 

■ »e B ( t 

t2 2 ) 

have been 

used. 




Multiplying both sides of Egs. (3.51) and (3.54) by ce n ( *l ) and Egs. 
(3.52) and (3.53) by se n ( l ), Integrating with respect to 1 froa 0 to 2z, and 
applying the orthogonality relations of the angular Mathieu functions leads to 


(3.63) 

*ln Ce n * J;* ^la^^a^B/n 

(3.64) 

®ln Se n * Jt' ®2a Ge ) c a a B,n 

(3.65) 

B jv. ^laCe.^a - v)iB ln Se n * * 


B T lh 2 /*2 2 Jv. *2a pe *a - W»Yih 2 /Y 2 2 J,' B 2 .Gek B '* n/B 

(3.66) 

w€ l A ln Ce n + ®la Se a , )'n,a * 


"Vlh 2 /»2 2 & W«ka'®n,a ♦ BT lh 2 /Y 2 2 B^Gek.^,. 


The priae over the suaaatlon sign is used to indicate that either odd or even 
values of a are used accordingly as to whether n is odd or even. 

«a,n/ Ba,n/ H*B, n and are glven by th * follo * ln 9 

(3.67) « 1/n * J^se B 4 (Dse n m dl / X*se n 2 (T) dl 

(3.68) 8 1/n « jTce B t (1)ce n ( , t) d7 /j£“ce n 2 (t) dl 

(3.69) 4y n « Jj 1 se 1 '(t)ce n (t) dl /£? l ce n 2 (l) dl 

»k,n (*)se n (1) dl /STi ie n 2 (1) dl 


(3.70) 
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Making use of Eqs. (3.(3) and (3.64), Eqs. (3. 65) and (3.(6) yields tvo 
sets of infinite homogeneous equations 

,, !k. * 2 * s, ' n * 5 .’. B2,t *- n • 0 

k * k ■»*.. • « 

vhere 

(3.72) g., n « -d - Y lh 2 /Y 2 2 )Fek B ( £ 0 ) 6 mfX » z , n 

(3.73) h Bfn ■ wH« >fn /8 ( Gek B ( E 0 )Se B '( l 0 )/Se B ( £ 0 ) 

- Gek B '(£ 0 ) Y lh 2 /Y 2 2 ) 

(3.74) s 1/n ■ vB >fn /8 I ^Fekat 6o> Ce m'( toJ/Ce,! £ 0 ) 

- C 0 Pek B ' (£ 0 ) 7 lh 2 /Y 2 2 ) 

(3.75) t B#n « (1 - Y lh 2 /Y 2 2 )Gek B ( E 0 ) ^ ' ««,r^r f n 


For a nontrivial solution/ the infinite determinant of Bq. (3.71) must 
vanish. The propagation constant 8 can then be determined from the roots of 
this infinite determinant. 

The infinite deteminant for odd values of m and n is 


(3.76) 


sll til s31 t31 - 
gll hll g31 h31 - 
sl3 tl3 s33 t33 - 
gl3 hl3 g33 h33 - 


0 


: 
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and for even values of a and n is 


sOO 

■ 20 

t20 

■ 40 - - 

s02 

s22 

t22 

■ 42 - - 

^02 

q22 

h22 

942 - - 

■ 04 

■ 24 

t24 

■ 44 - - 

t 

1 

1 

: - - 

1 

• 

■ 

• 

• 

: - - 



4. WEAKLY GUIDING APPROXIMATION 

The exact characteristic equations obtained in Chapter 3 are valid for 
an anisotropic elliptical fiber with any eccentricities. These equations are 
also applicable to the fiber with any refractive index differences between the 
core and cladding aaterial. However, for most of practical fibers, the 
difference in the refractive indices of the core and the cladding is typically 
very saall. The simplified characteristic equations can be obtained under 
this condition which is known as the weakly guiding approximation. 

Applying the weakly guiding approximation results in the following 

(4.1) y 2 2 - Y lh 2 ♦ Hwfycj (1 - c 0 / € l> * *lh 2 

(4.2) 1 - Y lh 2 /Y 2 2 * 0. 

4.1 EVEN MODES 

Applying Eqs.(4.1) and (4.2) into Eqs.(3.20) and (3.30) yields the 
equations 

(4.3) ce n ( 1 , t 2 2 ) « ce n ( t , T lh 2 ) 

<«.«) « jfce.l 1 , T lh J)ce„( 1 , 1 . ''lh 2 >« 

* Am,n 

where A, /n is the Kronecker delta which is zero when m t n and is unity when 

a ■ n. 

Substituting Equations (4.1) - (4.4) into Equations (3.35) - (3.38), the 
infinite determinants for even modes become 

(4-5) 1! ( 9a, a *m,m " h «,» 8 »,» ) “ 0 

or 

(4.6) 9a, a^, a “ h a,a 8 m,a * 0 

for a * 0, 1, 2, - - - . 
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By substituting Equations (3.35) - (3.38) Into Eq.(4.6), the following 
equation is obtained 

-d - Tih 2 / t 2 2 )( € l/€ 0 * ' r lh 2 / v 2 2)( l 5£! p »,n J/ n,B l &*8,n 

(4.7) ■ (Ce > , (6 0 )/Ce B (g o ) - (Yx^^JPek.* (£ 0 )/Fek B (£o) ) 

( (fix/fioJSe,' (E 0 )/Se B (£o) -(Y lh 2 /Y 2 2 )Gek B ‘ (£ 0 )/Gek B (Lo) ] . 
This is the siaplified characteristic equation for even nodes compared to the 
infinite deteralnants as given in Eqs.-(3.39) and (3.40). When the elliptical 
rod degenerates to a circular rod, the siaplified characteristic equation 
becoaes that of the anisotropic circular fiber. 

4.2 ODD MODES 

Applying Eqs.(4.1) and (4.2) in Bqs.(3.58) and (3.67) yields the 
equations 

(4.8) se n (^, Y 2 2 ) * se n (*l, Tj^) 

(4 .9) « B , n « iTse.Cl, Y lh 2 )se n <$, Y lh 2 )dg/ //%e. 2 (l, t lh 2 )d1. 

c A 

B,n 

where n is the Kronecker delta which is zero when a # n and is unity when 

a ■ n. 

Substituting Equations (4.1) - (4.2) and (4.8) - (4.9) into Equations 
(3.72) - (3.75), the infinite determinants for odd nodes become 

(4.10) 1 ( g B/ , t B/B - h B , B s B#B ) - 0 
or 

14.11) 9s,m^m,a “ h a,m s a,m * 0 

for a « 0, 1, 2, - - - . 

By substituting Equations (3.72) - (3.75) into £q.(4.11), the following 


equation is obtained 
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"Cl - Vlh 2 / T 2 2 >< £ l/ c o “ Y lh 2 / T 2 2 >< jfc P a,n *VaJi| a a,n Yn,a/ # a,a®B,a> 
(4.12) - lSe B '(l 0 )/8e«(io) " C »il» 2 /T2 2 > 0«k* * (f 0 >/Gek B (g 0 ) ] 

Kei/eoJCeaMEoJ/Ce^tSo) - < Y i h 2 /T 2 2 )?«*■' <£o>/ Fek a<io> 1 • 
This Is the slaplifled characteristic equation for odd aodes coapared to the 
Infinite deterainants given in Bqs.(3.76) and (3.77). When the elliptical rod 
degenerates to a circular rod, the slaplifled characteristic equation becoaes 
that of the anisotropic circular fiber. 



5. NUMERICAL RESULT8 FOR PROPAGATION CONSTANTS 


5.1 I8QTROPIC ELLIPTICAL FIBERS 

then the anisotropic parameter g in Eg. (3.1) is equal to unity, the 
siaplified characteristic equations in Bq.(4.7) and Eq.(4.12) becoae that of 
an isotropic elliptic guide. In Figures 2 through 5, the noraalized guide 
wavelength X/Xq for the isotropic elliptical fibers is plotted as a function 
of the noraalized cross-section area and noraalized aajor axis for €i/e 0 * 2.5 
and for the various values of £<,. These results are coapared with those given 
by ¥eh( 2 ) which are indicated by syabols; the results are in close agreeaent. 

For the eHE^i a ode, it can be seen in Figure 2 that the noraalized guide 
wavelength is alaost equal to unity for the very saall value of the cross- 
section area. This indicates that the geoaetry of the waveguide has no effect 
on the noraalized guide wavelength when the wavelength is such larger than the 
physical diaension of the core of fibers. For a fixed value of cross-section 
area, the noraalized guide wavelength is saaller for larger the value of 
This indicates that aore energy is carried inside of the circular core than 
the elliptical core. As the noraalized cross-section area becoaes larger, the 
difference in the noraalized guide wavelengths for varying | 0 becoaes saall 
again. This is since aost of the energy is carried inside of the core and the 
geoaetry of waveguide has no effect on the noraalized guide wavelength. 

However, as observed in Figure 3, the oHB^i aode is different froa the 
eKEii aode in that the difference in the noraalized guide wavelengths for 
varying | 0 is saaller than that of eHEji when the value of noraalized cross- 
section area is fixed. This saall difference is due to the fact that the 
electric lines are being coapressed such that the field density is aore 
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concentrated Inside the waveguide. For a fixed value of cross-section area, 
■ore energy is carried inside of the elliptical core than the circular core 
since the normalized guide wavelength is smaller for smaller the value of Iq. 

In Figures 4 and 5, the normalized guide wavelength is plotted against 
the nomalized major axis for various values of Eq and for &i/6 0 * 2 * 5 * In 
these figures, the difference in the normalized guide wavelengths for varying 
§0 is larger than those in Figures 2 and 3 since there is more binding 
dielectric material in a circular core than in a flatter rod (i.e. smaller to ) 
when the value of normalized major axis is fixed. 

5.2 ANISOTROPIC ELLIPTIC AL FIBERS 

In Figures 6 and 7, the normalized guide wavelength for an 

anisotropic elliptical fiber is plotted as a function of the normalized cross- 
section area for various values of anisotropy and for %/€o * an< ^ 

^ * 0.5. These figures indicate that the geometry of the waveguide and 
anisotropy of the core have no effect when the wavelength is much larger than 
the physical dimention of the core of fibers which indicate that most of the 
energy is carried outside of the core. For a fixed value of cross-section 
area, the normalized guide wavelength is smaller for larger the value of 
anisotropy. This condition indicates that the field intensity is more 
concentrated in the core, thus indicating that more energy is carried inside 
of the core. As the normalized cross-section area becomes larger, the 
difference in the normalized guide wavelengths for the varying anisotropy 
becomes smaller again. This indicates that the geometry and anisotropy of 
waveguide have a smaller effect on the normalized guide wavelength. 
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The noaallzed guide vavelength in Figures 0 and 9 is plotted against 
the nomalized aajor axis for the various values of anisotropy and for €i/€q * 
2.5 and Iq • 0.5. The effect of anisotropy on the normalized guide vavelength 
is sinilar to those in Figures 6 and 7. 









6. POWER CONSIDERATIONS 


The power along the z exie in the Bedim i of the fiber say be obtained 
by integrating the poynting vector over the surface area. 

Pi » 1/2 J f (E t X H t * ) • 2 ds 

- 1/2 - E^Hn*) L 2 d*d€ 


( 6 . 1 ) 

where 


L » ql( cosh 2§ - cos 2*1 )/21* 
and E 0 » 0 and t 2 * • . 


6.1 EVEN HOPES 

Substituting Equations (3.5) through (3.8) into Eq.(6.1) and integrating 
over the core area yields 

(6.2) Pcore ■ (1/2?!*) £°t 3**1 Jl^lm 2 * *Se B ' 2 ♦ S M Se B 2 > 

+ 0whJ: # B 1b 2 ( «Ce B ’ 2 4 C, B Ce B 2 >1 d£ 

4 ( 6 2 +v 2 HEi)/2Yi 4 (>( P ? Bi B Ai n T Bn l Ce B Se n f* 

♦ 0wH/2Y 1 4 X?.l BlaBlnCanlCea'^n ' C^n't^a - an)* 1 

4 0wEi/2Yi 4 Jc i ?' Ai B Xi n S Bn [Se B 'Se n - 8e»Se n * £*(!>, - b n ) 1 
Sinilary, the power carried in the cladding is 

(6.3) Pclad - (1/2y 2 4 )JJ*1 Bve 0 E, A 2b 2 ( «Gek B ’ 2 4 S BB *Gek B 2 } 

4 0wM i S # B 2b 2 ( *Pek B ' 2 4 C^Fek, 2 )) d! 

4 ( 0 2 4w 2 J1E o )/2Y 2 4 J^?^B 2B A 2n T Bn *( Fek B Gek n j^ 

4 0w)i/2Y 2 4 t ?^ B 2B B 2n C Bn *(Fek B 'Fek n 

- Fek B Fek n ' £ja B - a n ) _1 
4 0vE o /2Y 2 4 J? ( £1 A 2B A 2n S Bn *(Gek B , Gek n 

- Gek B Gek n ' ^(b, - b n ) _1 
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In Kqs. (6.2) and (6.3), a. and b, aze the characteristic values of the even 
and odd Kathleu functions of order a, respectively. The prise over the 
sanation sign is used to Indicate that ■ » n is excluded. Also, Yj 2 ■ 
4Yih 2 /q 2 and Y 2 2 ■ 8 Y 2 2 /q 2 are used. 

The following abbreviations have been used. 


(6.4) 

Can 

■ J? ce B ' (1)ce n ' CJ) d2 

(6.5) 

*mn 

■ £^e B '(t)se n '(1) dt 

( 6 . 6 ) 

T an 

* ifce B ' (?)se n (1) dt 

(6.7) 

C.n* 

* iTce B ,, (7)ce n *'( , l) dl 

( 6 . 8 ) 

San* 

■ C^e/'CtJse^Mt) da 

(6.9) 

Tan* 

- jf ce B # *(t)se n *(a) da . 


The power distribution characteristics for the eHEn node is given in 
Figure 10. The fractional power carried by the core and cladding is plotted 
against the noraallzed najor axis for the various values of anisotropy and for 
C ^/£ 0 ■ 2.5 and t 0 » 1 . 0 . Most of the power is carried in the cladding near 
the cut-off and in the core far froa the cut-off. For a fixed value of the 
nornalized aajor axis, the nore energy is concentrated inside of the core for 
larger the value of anisotropy and far fzon the cut-off. 

6.2 OOP HOPES 

Substituting Equations (3.43) through (3.46) into Eq.(6.1) and 
Integrating over the core area yields 

(6.10) Pcore « (1/2Y 1 4 )^? 0v E A 1b 2 ( xCe B ' 2 ♦ C^Ce * 2 ) 

♦ PvH E Bi a 2 { xSe B ' 2 ♦ S BB Se B 2 )] d l 
+ ( • 2 *w*|ie 1 )/2Y 1 « J5 5 # B lB A ln T Bn l Se B Ce n J° 

♦ 6wp/2Y 1 4 | ? e e; B lB B ln S Bn (Se B ‘Se n - Se B Se n * {*(!>, - b n r 2 

♦ 0wc 1 /2Y 1 4 JS 1* li.AmCsnlCea'On " Ce B Ce n '(°(a B - a n ) _1 
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Simllary, the power carried In the cladding is 

(6.11) Pclad ■ (1/2Y 2 4 > evc 0 j? rfr i 2 « 2 < «Fek B '2 ♦ C./Fek* 2 } 

♦ 0v)l E B 2b 2 * «Gek B ' 2 ♦ 8 M 4 Gek B 2 )) d£ 

♦ ( e 2 *w 2 MC 0 )/^ 2 4 j 5 ,£ e B 2B * 2n T Bn t ( Gek B Fek n £ 

♦ evM/2T 2 4S # e; B 2l B 2n S >n 4 lGek > *Gek n 

- Gek a Gek n ' £(b a - b n )'l 
+ Bve 0 /2Y 2 4 ? ( e; k 2a i 2n C an t (Pek B 'Fek n 

- Pek B Fek n '£(a B - a n )-l 

In Sqs.(6.10) and (6.11), a B and b B are the characteristic values of the even 
and odd Hathleu functions of order a, . respectively. The prlae over the 
suaaation sign is used to indicate that a » n is excluded. Also, y 2 2 ■ 
4Yijj 2 /q 2 and Y 2 2 « 4Y 2 2 /q 2 are used. 

The following abbreviations have been used, 

(6.12) « Jjfce B ' (*)ce n * (a) da 

(6.13) S Bn * jfse B ' (1)se n ' (1) da 

(6.14) T Bn * Xfce B (a)se n '(a) da 

(6.15) C.H* « J?’ce B , ’(t)ce n t '(a) da 

(6.16) S Bn * - j^'se B *'(a)se n t '(a) da 

(6.17) T an * - J^ce B *(a)se n t ' (a) da. 

The power distribution characteristics for the oHEu aode is given in 
Figure 11. The fractional power carried by the core and cladding is plotted 
against the normalized aajor axis for the various values of anisotropy and for 
%/£<> “ 2 * 5 * nd *o " 1*0* Most of the power is carried in the cladding near 


the cut-off and in the core far fros cut-off. However, the difference in the 
power distribution for oHB 22 for the varying anisotropy is ssaller than that 
of eHEn for a fixed value of normalized sajor axis. 





7. DESCRIPTION OF COMPUTER PROGRAMS 

In this chapter, the computer programs for the anisotropic elliptical 
fiber vill be considered. These programs are written in the language of 
FORTRAN IV and in order to present the complete computer programs the 
subroutines to calculate the Mathieu and modified Mathieu functions developed 
by Rengarajan and Levis(49] vill be included. The theory and notations used 
in the computer programs are the same as those employed by McLachlan[48] . 

The normalized propagation constant and power distribution 
characteristics as a function of the normalized cross section area or major 
axis for the given value of and anisotropy have been determinded by 
utilizing these programs. These computer programs consist of a main program 
and user called subroutines: CHAVAL and POWER. These subroutines CHAVAL and 
POWER call nine subroutines in order to compute the Mathieu and modified 
Mathieu functions. 

In subroutine CHAVAL, an initial guess for the given mode is chosen and 
used to evaluate either Eg. (4.7) or Eg. (4.12). Next, Muller's method is used 
iteratively to determine the normalized wavelength that vill minimize the 
function; an error criterion has been used to terminate the Iteration. In 
subroutine POWER, the power distribution characteristics are calculated using 
the normalized wavelength obtained in subroutine CHAVAL. The algorithm was 
run on an CYBER 990 using only a single processor. 

The seguence of the called subroutines is illustrated in the following 
flow chart. 
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8. CONCLUSION 


The exact characteristic equation for anisotropic elliptical optical 
fibers is obtained for the odd and even hybrid modes in teras of infinite 
determinants employing Hathieu and modified Mathieu functions. The exact 
characteristic equation Is applicable to elliptical fibers vith any 
elllpticity. A simplified characteristic equation can then be obtained by 
applying the veakly guiding approximation such that the difference in the 
refractive Indices of the core and the cladding is small. Under this 
approximation, it can be shovn that significant simplification can be 
achieved. 

The simplified characteristic equation is used to compute the normalized 
wavelength for an anisotropic elliptical fiber. When the anisotropy parameter 
Is equal to unity, the characteristic equation becomes that of Isotropic 
fiber. The results are compared to the previous research and they are in 
close agreement. Por a fixed value of the normalized cross-section area or 
major axis, the normalized wavelength >/A 0 is small for larger the value of 
anisotropy. This condition indicates that more energy is carried inside of 
the fiber. However, the geometry and anisotropy of fiber have a smaller 
effect when the normalized cross-section area or major axis is very small or 
very large. 

An exact solution for the wave equation can not be determined when the 
thermoelastic stress causes a transverse anisotropy over the core of fibers. 
One possibility is that the propagation characteristics in the biaxial 
anisotropic fibers could be obtained by applying the numerical or 
approximation techniques given in Chapter 1 and this could be a subject for 
further study. 
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APPENDIX 

The following is a listing of the computer programs, MAIN 
POWER, CHVAL2, EXPAND, ANGMFC, FACTOR, 8TORE, CERAD, 8ERAD, FBRAD 
written in FORTRAN IV language. 


CHAVAL, 
AND GERAD 
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C 

C THIS PROGRAM CALULATES THE NOMALIZED WAVELENGTH AND 

C POWER DISTRIBUTION FOR ELLIPTICAL FIBER AS FUNCTION OF 

C NORMALIZED CROSS-SECTION AREA OR MAJOR AXIS* 

C 

PROGRAM MAIN 

IMPLICIT OOUBLE PRECISION <A-H,0-Z) 

DIMENSION RES ( 56 1 f X(56) 

OPEN <6»FILE*‘0UTPUT* I 
C 

C NEVOO * * 1 FOR ODD MODE * 

C « 2 FOR EVEN MODE* 

C XASc * « l FOR NORMALIZED CROSS-SECTIONAL AREA, 

C ■ 2 FOR NORMALIZED MAJOR AXIS* 

C ETA : INDEPENDENT VARIABLE IN MATHIEU FUNCTION 

C PHI : INDEPENDENT VARIABLE IN MODIFIED MATHIEU FUNCTION 

C MODE t WAVE MODE NUMBER 

C P : EP1/EPO 

C G : ANISOTROPY EPZ/EPX 

C 

c 

NEVODM 
K AS E *2 
MODEM 
P *2 • 5D0 
GM.5D0 
PHI *1 •CO- 

GO TO 111,12), NEVOD 

11 WRITE(6,131) MOOE 
GO TO 13 

12 WRITE(6fl55) MODE 

13 WR1TE(6,132) P 

WRITE(6,1D3) S 
WRITE(6,104) PHI 

C 

C CALCULATE THE NORMALIZED WAVELENGTH 

C 

CALL CHARVAL C NEVOD , P , G, MOOE »PH I , X AS E , R SS I 
C 

DO 14 1*1,56 

14 xm*REsm 

c 

C CALCULATE POWER DISTRIBUTION FOR THE GIVEN MODE 

C 

CALL POWER (NEVOD, P,G, MODE, PHI ,KASE,X) 

C 

121 FORMAT ( ‘ 1 ‘ » ‘THIS IS RESULT FOR ODD MODE, M *•, 12) 

122 FORMAT (IX, ‘RATIO OF CORE AND CLADDING PERMITTIVITY *‘, 

* 012*5) 

123 F0RMATI1X, ‘ANISOTROPY *», D12.5) 

i:4 F0RMATC1X, ‘VARIABLE IN MATHIEU FUNCTION *‘, D12.5) 
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1C5 FORMAT ( * l ' * • THIS IS THE RESULT FOR EVEN MODE* M«*. 121 
STOP 
E ND 

SUBROUTINE CHARVAL ( NE VOD* P * G» MODE »PHI*KA$E»RESI 


PURPOSE * CALULATE THE NORMALIZED WAVELENGTH FOR 
ELLIPTICAL GUIDE. 

INPUT s NSVOO -(INTEGER) SPECIFIES ■ 1 FOR °jj EN M jjjjp E 

P -(DOUBLE PRECISION) IS THE RATIO BETWEEN CORE 
AND CLADDING PERMEABILITY. 

G -(DOUBLE PRECISION) IS THE ANISOTROPY* EZ/EX. 
MODE -(INTEGER) IS THE MODE OF CHARACTERISTIC 
EQUATION. 

PH! -(DOUBLE PRECISION) IS INDEPENDENT VARIABLE 
IN MODIFIED MATHIEU FUNCTION. 

XASE -(INTEGER) * 1 FOR NORMALIZED CROSS-SECTION 

AREA* 

« 2 FOR NORMALIZED MAJOR AXIS. 
OUTPUT * RES -(DOUBLE PRECISION) CONTAINS THE NORMALIZED 

WAVELENGTH. 


C 

c 

c 

c 

c 

c 

c 


21: 


c 

c 

c 


IMPLICIT DOUBLE PRECISION (A-H»0-Z) 

DIMENSION CHVK23), CHV2I23), AB(25), QV « ♦ 1 • ET A 1 21 » • 

* SEl(21)»SElD(21)»SE0(21)*CEi(21)*CElD(21)» 

* CEC(2l)»CS0i(21 )* SEC 1(21 ) . CEDSEC 21 ) » SEDCE 1 2 1 ) * 

> SE 1 S Q ( 2 1 )*CE1SQ(21 ) ,CEGM( 21 ) » SEOM( 21 ) *CE1M( 21 ) * 

» SE1M( 21 ) »RES( 56) 

DATA P 1/3. 1*1592653539793 DC/ 


c t A : INDEPENDENT VARIABLE IN MATHIEU FUNCTION 

XI s INDEPENDENT VARIABLE IN MODIFIED MATHIEU FUNCTION 
I ORDER : WAVE MODE NUMBER 
p : EP1/EPC 

G s ANISOTROPY EPZ/EPX 


M*2C 

last»*o 

I F ( KASE « EQ« 2 ) LAST»56 
X I »PHI 

DT*0TANH(XI) 

DCOS2*DCOSH(XI)*DCOSH(XI ) 

X2=PI*DT*DCOS2 
WRITE(6.210) X2 
FORMATdX* 012.5) 

SUBDIVIDE ETA FOR INTEGRATION. ALFA, BETA, GAMMA AND ANU. 

H»2.000*PI/20.CDC 
£TA( 1 )-0.0D0 
00 11 1 = 2*21 
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1 STAC I )*ETA( 1 )♦< I-l |*h 

CSA - NORMAL I ZED CROSS SECTION AREA OR MAJOR AXIS 

CSA«0.000 
X I NC*2 . 50-2 
V ARX*0. 99999999 90C 
DO 70 K*l * LAST 
CSA*C$ A+X I NC 
GO TO (12.13), KAS E 

12 WRITE(6,201) CSA 

201 FORMATUX, ‘NORMALIZED CROSS SECTION AREA « ‘.D12.5I 

X1«CPI*pi#cSAI/(4.0DC*DT*DC0S2> 

GO TO 14 

13 HR l TE ( 6 » 2C 2 ) CSA 

2-2 FORMAT! 1 X» ‘NORMAL I ZED MAJOR AXIS « ‘*012.51 

X1»<PI*PI*CSA**2)/C4.0DC*DCOS2I 

c 

^ FIRST GUESS OF VAR * LANDA/LANDAO 

c 

1A IFCK.EO.ll VAR» VAR X 

NC* 1 
NS*Q 

1C VAR2*1 .C DO/ ( VARPVAR ) 

C 

c EVALUATE 3 s INDEPENDANT VARIABLE 

C GAMMAE * QVCll, GAMMAH * Q V ( 2 ) * GAMMA2 * 0VC3) 

OVC 1 )*X1*( P*G-G*VAR2 ) 

QV!2)*X1*(P-VAR2» 

0VC3I*X1*C1.0D0-VAR2) 

QV(4)*X1*( 1.000- VAR 2 I 

c 

c CALCULATE MATHIEU AND MODIFIED MATHIEU FUNCTIONS 

DO 53 KQ* 1*4 

GO TO (15*16)* NEVOD 

15 I EV0D*1 

IF(M0D(K9*2).EQ.1 ) IEV0D*2 
I ORDER*MODt 

I F ( MQDC KQ* 2 ) .EO.l ) IORDER*MODE 
Q*9V( KQ) 

CALL CHVAL2 CM* Q*CHVl *CHV2 » J ) 

CV*CHV1 C IOROER ) 

I FCMODC KQ* 2 ) .EO. 1 ) CV*CHV2 ( IORDERM ) 

GO TO 17 

16 1EV00*1 

IF(MOO(KQ*2).EQ.O) IEV0D*2 
I 0R9ER*M0DE 
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c 

c 

c 

17 


21 

41 

23 

42 

22 

43 

24 

44 
C 
C 

c 

45 
C 
C 

c 

c 

c 

c 

51 

53 

52 

54 


IFCMODCKQ,2)«EQ«C) !ORDER*MODE 
Q-QV(KQ) „ , % 

CALL CHVAL2 C M* Q tCHVl »CHV2 » J ) 

SJSSoUqSEU.o* 

OBTAIN EXPANDING COEFF1CI ENT . ABXX 

CALL EX'PANDCQ*IEVOD*IORDER*CV*3»ABtN> 

CALCULATE MATHIEU FUNCTIONS AND DERIVATIVES, 

OROER ■ MODE 

KO c O*KQ 

IFCMODCNEVOD,2).EQ.l) KQEO-KQ+4 
GO TOC21. 22. 23. 24, 22.21,24, 23). *Q E0 

SEIII >*»NSHFCIO 1 ievO0.IORDE».ET»(l 

SE10m*ANGMFCC0. !EVOD#IORD€R. ETAm.lt AS. N» 

GO TO 45 

S Eot 1 1 *ANGHFC( QtlEVOOtlOROERtETAI ll*0»A8*N) 

GO TO 45 

CElCIliANGMFCCO.IEVOD.IORDER.ETAf Ij.O.AB.NJ 

CE1DC I l*ANGMFCCQ. IEVOD. IORDER.ETAC I ) .l.AB.N) 

GO TO 45 

CEO( ni*NGHFCI#»ISVOO.IO*DER.ST»l 1 1.0. AS. N) 

NORMAL ! ZAT ION FACTOR for HOOIFIEO MATHIEU FUNCTION 

CALL FACTOR ( IEVOD. I ORDER » 0. AB. N.PS I 

COMPUTE AND STORE THE VALUES OF BESSEL FUNCTIONS 

CALL STORECQ.Xl.N) 

CALCULATE MODIFIED MATHIEU FUNCTIONS 

GO TO C51.52.53. 54, 52.51. 54.53). KOEO 
SE«SERAOCQ» I ORDER. 0 » PS. AB.N) 

SED*SERADCQ» I0R05Rtl.PS. AB.N I 

GO TO 50 „„ 4a M . 

ge*geradco.iorder.o.ps.ab.n) 

GED«GERADCQ. IOROERtl.PS. AB.N) 

C E*CER ADC Q»1 OR DER.O. PS, AB.N) 

CED«CERADCQ»I0RDER.1.PS.A8,N) 

FE«FERADCQ, 1 OR DER.O. PS tAB.N) 


nonnu» 
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FEO«FERAO!Q» I ORDER . 1 » PS » A9 »N 1 
CONTINUE 

CALCULATE M TH TERM! « MODE I OF 
ALFA* BETA* GAMMA AND NU. 

WRITE! 6.1071 S E , S EO . G5 . GE 0 , C E . CEO. FE »FE 0 
00 56 I «1 » 2 1 
CEOMC I I *CEO( I ) 

SEOMC I 1-SEO! I ) 

CE1M! 1 I*CE1! I 1 
sEiwm-SEim 
CEOim*ceom*CEiui 
SE01( I l*SEO( I )*SE1! I » 

SE1SQ! I )>SE1U)*SE1( 1 I 
CE1SQ! II-CE1! I1*CE1! I 1 
CEDSE! I 1 *CEID! I I*$E1! I I 

56 SEDCE! I 1*SE1D!I 1*CE1! 1 1 
$1«SIHPSN!CE01»20*H> 

S2«SIMPSN!SE0l,20»H) 

S3«SIMP$N!CED$E,2C»H) 

S4«SIMPSN!$EDCE.20»H> 

S5«S1MPSN!SE1S0.20»H) 

S6*SIMPSN!CE1SQ»2C .H) 

S 5M*S 5 
S6M*S6 

GO TO 157.581. NEVOD 

57 ALFAM*S2/S 5 
BET AM*S1/S 6 
GAMMAM«$4/S6 
ANUM*S3/S5 
GO TO 59 

58 ALFAM*S1/S6 
BET AM* S2/S 5 
GAMMAM*S3/S5 
ANUM«S4/S6 

59 XMSQO«ALFAH*BETAM 
RHC*QV!2I/QV!3) 

WRITE! 6. IC3) AL F A M , B 6TAM , GAMMAM . ANUM , RHC 
1 C 3 FORMAT! IX. 5012.5) 

C 

C CALCULATE MATHIEU FUNCTION INTEGRALS 

C OROER ■ N. N*2 

C 

ALFA«O.CDO 

BETA«0.000 

GAMMA *0.0 DC 

ANU-0.000 

XMSQN1«0.000 

XMSON2«O.OOC 

00 90 I M* 1 * 4 
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I 0RDcR*2* I M-l 

IFIM00«M0DE*2).EQ.C > I0RDER«2*IM-2 

IFIIOROER.NE.MODEI GO TO 61 

ALFA*ALFAH 

8ETA*BETAM 

GAMMA*GAMMAM 

ANU*ANUM 

XMSQN1*XMSQN1+3ETA*ANU 
XHSQN2*XMSQN2+ALF A*GAMM A 
XMSQ«f XMSQN1*XMSQN2) /XHSOO 
GO TO 92 

61 DO 80 KQ*1 * 4 

GO TO (62*631* NEVOO 

62 I EVOD* 1 

IF(M0D(KQt2I.EQ.l I IEV0D«2 
Q*QV(KQ» 

CALL CHVAL2(H,Q,CHV1 ,CHV2,J) 

CV-CHVl ( IOROERI 

IF(MOD(KQ*2».E0.1 » CV*CHV2 I I 0RDER*1 » 

GO TO 64 

63 I EVQD*1 

IF(MOD<KQ,2I.EQ.OI IEV0D*2 
Q*QV(KQI 

CALL CHVAL2(M»Q»CHV1*CHV2*J) 

CV*CHV1 C IOROER) 

1FCM0D(KQ*2I.EQ.0» CV«CHV2( I ORDERED 
C 

C OBTAIN EXPANDING COEFFICIENTt ABXX 

C 

64 CALL EXPANDIQ* IEVOD* I0RDER»CV,3»A8»N) 

C 

C CALCULATE MATHIEU FUNCTIONS AND DERIVATIVES 

C 

KQEO-XQ 

IF(M0D(NEV0D*2I .EQ.l I KOEO*KQ+4 
GO T0(71*72»73»74*72»71,74.73l, KQEO 

71 DO 81 1*1 t 2 1 

SE1C I l*ANGMFCf QtlEVODt I ORDER *ETA ( 1 l,C*AB*N) 
81 SE1DI I I * ANGHFC ( 0 ♦ IEVOD* IOROER . ETACII • 1 »AB *NI 
GO TO 8C 

73 DO 82 1*1 » 2 1 

92 SEOC I l*ANGMFC(0» I EVOD* IOROER? ETAC I > *C * AB*N) 

GO TO 8C 

72 DO 83 1*1*21 

CE1 ( I )*ANGHFC(Q» 1EV0D* I ORDER * ET A ( I ) * 0 * AB * N ) 

8 3 CE1D< I ) * ANGHFC ( 0 * I E VOO* IORDE R » ET A * I I * 1 * AB » N ) 
GO TO 60 

74 00 84 1*1*21 

84 CEO (I I* ANGHFC (0* IEVOD *1 OR DER*ETACI)»0*A3*NI 

80 CONTINUE 

C 
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C CALCULATE SUM (BETA * NU> AND SUM (ALFA * GAMMA) 

C 

DO 86 1*1 *21 

CEOHl l«CEOM(I )*CE1( I ) 

SE01( I l*SEOM( I )*SE1(I) 
seisq( n-SEK n*sEim 
CElSQCI)*CEim*CEl( i » 

CEOSE( 1 )*CE1D( I )*SE1MC I ) 

86 SEDCE( I )*SE1D( I )*CE1M( I ) 

S1*S1MPSN(CE01»20»H) 

S2*$IMPSN(SS01*20,H> 

S3«SIMPSN(CEDSE»20»H) 

S4»SIMPSN($EDCE*20*H) 

S5*SIMPSN(SEIS0»20»H) 

S6*SIMPSN(CE1SQ*20*H) 

C 

GO TO (87,88)* NEV3D 

87 ALFA*$2/S5 
BETA*S1/S6 
GAMMA*S«/S6M 
ANU*S3/S5M 
GO TO 89 

83 ALFA*S1/S6 

3ETA«S2/S5 
GANMA*S3/S5M 
ANU*S4/S6M 

89 XMSQN1«XMSQN1*BETA*ANU 

XMSQN2*XMSQN2«>ALFA*GAHMA 
XMSQ*( XMSQN1*XMSQN2 )/XMSOD 

92 WRITE(6,103) AL F A , BET A , GAMH A , ANU , XMS Q 

9C CONTINUE 

C 

C EVALUATE CHARACTERISTIC EQUATION 

C 

Y1*-(XMSQ*(1.0D0-RHC)**2) 

Y2*VAR*VAR 

GO TO ( 93 . 94 ) » NEVOD 

93 Y3*(SED/SE)-(RHC*GED/GE) 

Y4«(P*CED/CE)-(RHC*FED/FE ) 

HR I TE ( 6 , 10 1 ) Yl, Y2* Y3» V4 
GO TO 95 

94 Y3»(CED/CE)-(RHC*FED/FE) 

Y4*(P*SED/SE)-(RHC*GE0/GE ) 

WRITE (6, 101) Yl* Y2 » Y3, Y4 

95 YX*Y2*Y3«Y4 
YZ*YX/Y1 
Y*YZ-1.0D0 

WRITE (6*203) YX, Yl* YZ* Y, VAR 

213 FORMAT( 1X*5D12. 5) 

C 

C DESIDE ON TOLERANCES 
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C 

IFINC.NE.l.AND.DABSI Y).LE.2.00-3) GO TO 39 
IFINC.EQ.l) GO TO 32 
1 FI NC« EQ. 2 ) GO TO 34 
IFINS.EQ.O) GO TO 34 
IF! Y*YS 1 1 36,36,31 

31 YS1«Y 
VARS1«VAR 

VAR«C VARS1+VARS2)/2.0DC 
NS«NS*1 

IF4NS.IE.20) GO TO 10 
GO TO 47 

C 

C 1 ST CALCULATION OF Y, DECREMENT VAR BY O.C1 

C 

32 YS1«Y 
VARS 1 *VAR 

33 VAR«VAR-1 .00-2 
NC*NC* 1 

YMI N«Y 
VARMIN*VAR 

I F( NC.LE.2C ) GO TO 10 
GO TO 48 

34 IFCY4YS1) 36,36,35 

35 IFIDA3SCYI-0ABS IYS1 1 1 32»33»33 

36 YS2«Y 

V AR S2 * V AR 

VAR*C VARS1+VARS2I/2.C00 
NS*NS+1 
YM I N*Y 
VARM I N*VAR 

I FI NS.LE.20 ) GO TO 10 

47 WRITEI6,1C6) 

WR I TE ( 6 , 10 8 ) YMIN, VARM I N 
GO TO 39 

48 WR I TE I 6 » 10 2 ) 

V AR«VARX 

WR I TE I 6, 10 3 ) YMIN, VARMIN 
39 RESI K) *VAR 

7G CONTINUE 

RETURN 

101 F0RMAT(1X,4D12.5I 

102 FORMATC IX, 'ERROR s RESULT HAS SAME SIGN FOR 10 TRIES') 

1C6 FORMATI IX, 'ERROR * 1C TRY FAILED TO OBTAIN RESOLUTION') 

1C7 FORMATI IX, 3012. 5) 

103 FORMATUX, 'OBTAINED RESOLUTION D12.5, 

C 'MIN CALCULATED NORMALIZED WAVELENGTH *', 012.5) 

END 

SUBROUTINE POWER I NE VOD, P , G, MODE , BOUND , KASE , A ) 

C 

C PURPOSE : CALULATE POWER DISTRIBUTION ON ELLIPTICAL 
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GUIDE. 

INPUT * NEVOD -(INTEGER) SPECIFIES » 1 FOR ODD MOOE 

* 2 FOR EVEN MOOE 

P -(DOUBLE PREC IS I ON) IS THE RATIO BETWEEN CORE 
AND CLAODING PERMEABILITY. 

G -(DOUBLE PRECISION) IS THE ANISOTROPY* EZ/EX. 
MODE -(INTEGER) IS THE MOOE OF CHARACTERISTIC 
EQUATION. 

PHI -(DOUBLE PRECISION) IS INDEPENDENT VARIABLE 
IN MODIFIED MATHIEU FUNCTION. 

A -(DOUBLE PRECISION) IS THE NORMALIZED 
WAVELENGTH. 

KASE -(INTEGER) * 1 FOR NORMALIZED CROSS-SECTION 

AREA* 

* 2 FOR NORMALIZED MAJOR AXIS. 
OUTPUT : RES -(DOUBLE PRECISION) CONTAINS THE RATIO OF 
POWER DISTRIBUTION. 


C 

C 

C 

C 

C 

C 

C 

C 


IMPLICIT OOUBLE PRECISION (A-H.O-ZI 
DIMENSION CHVK23I* CHV2I23)* A3(25l* QV(4), A(56l 
ETA(21I*PHI(41I*SE1(21I»SE10(21I*SS0(21) 
SE0D(21I*CE1(21)*CE1D(21 I .CEO! 21 ) .CECOI2 
SE ( 2 1 ) * S EDI 2 1 ) » CE( 21 ) » CED(21)» 

FEU1)» FEDUll* GEU1)* GEO(41)» 

S 1 ( 21 ) * S2 ( 21 ) * S 3 ( 2 1 ) * S4(21)» S5(2D* 
ST ( 2 1 ) » S8 ( 21 ) * S9I21I* S10(21). SI 1 ( 21 ) 
S 1 2 ( 2 1 I * S 1 3 ( 2 1 ) ♦ S14I21I* 

S21 ( 21 ) • S2ZC21). S23 ( 2 1 ) 

DATA PI/3.1415926535B97Q3D0/ 
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ETA : INDEPENDENT VARIABLE IN MATHIEU FUNCTION 
XI S INDEPENDENT VARIABLE IN MOOIFIED MATHIEU FUNCTION 
I ORDER : WAVE MODE NUMBER 
P * EP1/EPO 

G : ANISOTROPY EPZ/EPX 


M«2G 

EP*ll.OD- 91/(36.0 DO *PI) 

XNU«4.0 DC* PI *1.00-7 
CONS*DSQRT(EP/XNU) 

LAST*40 

1FIKASE.EQ.2) LAST*56 
PHI 0*0.000 
PHI1-B0UN0 
PHI2*5.CDC*3CUND 
OT«DTANH( PHI 1 ) 

DC0S2«DC0SH(PH1 1 > *DCOSH< PHI 1 ) 

X2»PI*DT*DCOS2 

WRITE (6*210) X2 

SUBDIVIDE ETA AND PHI FOR INTEGRATION. 
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C 

H1*2*C DO 4P 1/20*000 
STAC 1 1 ■0*000 
00 11 1 * 2 » 21 

11 STAC I l*ETAC 1-tl+Hl 
H2*30UNO/20.000 
PHI ( 1 1 *0*000 

00 12 1*2,21 

12 PHICn*PHl( I-1I+H2 
H3*CPHI2“PH! 11/20*000 
DO 13 I *22 » 41 

13 PHIIII*PHH1-1I*H3 
C 

c CSA - NORMALIZED CROSS SECTION AREA OR MAJOR AXIS 

C 

CSA*0.000 
XINC*2.5D-2 
00 70 K* 1 » L AST 
CSA*CSA»XINC 

IF(A(KI*EQ*l«ODO) CO TO 81- 
GO TO 114*151* KASE 

14 WRITE(6*201 1 CSA 

2C1 FORMAT ( 1X»* NORMALIZED CROSS SECTION AREA * *.012.51 

C 

Xl*tPI*PI*CSA!/<4.0DC*DT*DC0S2 1 
GO TO 16 

15 MRITEC6.202I CSA 

202 FORMATC IX, 'NORMALIZED MAJOR AXIS * *. 012.51 

C 

X1*<PI*PI*CSA**21/C4.0D:*DCCS2 1 

c 

c CALCULATE CONSTANTS. 

c 

16 VAR2*1.COO/(ACK1*A<K»I 
C 

C EVALUATE Q * INOEPENDANT VARIABLE 

C GAMMAE * QVCil, GAMMAH * QV ( 2 1 * GAMMA2 * OV ( 3 1 

C 

QVt 1 I*X1*CP*G-G*VAR2 I 
QVC2 1*X1*IP-VAR21 
QVC31«X1*C1.000-VAR21 
QVC41*X1*I 1.0DD-VAR2 1 
C 

C 4*C0NS / A C K I 
C1*P*C4 

C2*l*0D0/fA(K I * CONS I 

C3«P*VAR2 

C5*1.000»VAR2 

C6* ( QVC 2 1«QVC2 )I/IQVC3>*0VC3il 
C CALCULATE MATHIEU ANO MODIFIED MATHIEU FUNCTIONS 
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C 

00 50 KQ*l»4 

GO TO 117,181, NEVOD 

17 IEVOD-1 

1 F I MODI KQ « 2 ) • EQ* 1 ) IEV0D*2 
IORDER*MODE 

IF(M0DCKQ,2I.EQ.l I IORDER*MODE 
Q*OV IRQ) 

HR I TE ( 6 ,30 1 ) Xl »0 

CALL CHVAL2<M,Q,CHV1 ,CHV2, JI 

CV-CHV1I I0RDER1 

IFfM0DIKQ,2).EQ.l I CV*CHV2 1 1 0RDEIU1 1 
GO TO 19 

18 I EVOD*l 

IF(M0D(KQ,2).E0.C > IEV0D«2 
I ORDER»MOOE 

IFIMOD(KQ,2).EQ.O) I ORDER*MODE 
Q*QVUQ> 

HR I TE ( 6 * 3C 1 ) XI t Q 

CALL CHVAL2IM,Q,CHV1,CHV2»J> 

CV*CHV1 ( I ORDER ) 

IF(MODIKQ»2>.EQ.OI CV-CHV2 f I CRDER+l I 
C 

C OBTAIN EXPANDING COEFFICIENT, ABXX 

C 

19 CALL EXPANOIQ, I EVOD, I OR DER , C V , 3 . AB , N ) 

C 

C CALCULATE MATH I EU FUNCTIONS AND DERIVATIVES, 

C ORDER * MOOE 

C 

KQEO*KO 

IF(M00(NEV0D,2J.E0.i> KQE0*KQ«-4 
GO T0< 21, 22,23, 24, 22, 21, 24, 23), KQEO 

21 DO 41 1*1,21 

SE1(!)*ANGMFC(Q,I£VOD,IORDER,ETA( I ),C,AB,N) 

41 SE1DCI I *ANGMFCI Q»IEV0D,10RDER»ETA(IJ,1,AB,NI 
GO TO 45 

23 DO 42 1*1,21 

SEOOI I l*ANGMFCf 0,1 EVOD, I OR DER , ET A (I) , 1 , AB ,N » 

42 SEOC I »*ANGMFC(Q,1EV0D.I0RDER,ETA( I ),9,AB,NI 
GO TO 45 

22 00 43 1*1,21 

CE1( I l*ANGMFC(0, IcVOO, I OR DER, ETA I I 1 ,0 , A3,N) 

43 CE1DI I»*ANGMFC( Q, IEVOO, IORDER , ETA( I I ,1 ,A6,N) 

GO TO 45 

24 DO 44 1*1,21 

CEOD( 1 ) *ANGMFC ( 0, I E VOO, I OR DER , ET A (II , 1 , AB, N I 

44 CEOC 1 l*ANGMFC 10, IEVOO, IORDER »ET A ( I ) , D , AB, N ) 

C 

C NORMALIZATION FACTOR FOR MODIFIED MATHI EU FUNCTION 

C 
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45 CALL FACTOR! IEVOD* I ORDER , Q, AB , N,PS I 

C 

C CALCULATE MODIFIED MATHIEU FUNCTIONS 

C 

GO TO f 31, 32, 33. 34, 32. 31, 34, 331, KQEO 

31 DO 51 1*1,21 

CALL STORE l QtPH I< I l*N» 

ScCn-SERADCQ, IORDER.O,PS,AB,N> 

51 SED( I l«$ERAD(Q, IORDERtl »PS ,AB,N) 

GO TO 50 

33 DO 52 1*21,41 

CALL STOREIQ.PHlf I ) ,N I 

GEC I1*GERADIQ,10RDER,0,PS,A3,NI 

52 GED( I ) »GER ADI 0* 10RDER , 1 » PS , AB,N> 

GO TO 5C 

32 DO 53 1*1 *21 

CALL STORECQ*PHHt»*NI 
C E( I l*CERADI 0* !ORDER»0*P$ * AB.N) 

53 CED( I )*CERAD( Q* IORDERtl , PS, AB,N» 

GO TO 50 

34 DO 54 1*21*41 

CALL STOREIQ*PHim »N) 

F E I I I*FERADIQ,IORDER,0,PS,A3,N1 

54 FED! I l*FERADIQ, IORDERtl .PS.AB.NI 

50 CONTINUE 

C 

C CALCULATE INTEGRAND 

C 

DO 56 1*1*21 

S21I I »*SEO! !»*SEl( I ) 

S221 I I *SE 1 < I )*SE1 I I ) 

S23C I I *CE10( I I*SE1 I I ) 

SKI )*SE1D< I I*SE1D(I) 

S 2 ( I l*CE10( I >*CE1DI I I 
S3C I )*CclD( I I*SE1 1 1 I 
S 4 ( I )*SE0D( I l*SEODC I I 
S 5 C I )*CEOD( I )*C EO 0 1 I » 

56 S6( I l*CEODC I I *S EO ! I ) 

C 

C 

00 57 1*1,21 

S7C I ) *SED( I l*SEDI I ) 

SSI I >*SEU l«SE( I ) 

S 9 1 I l*CED( I )«CED( I I 

57 S10I 1I*CE( n»CE(II 
C 

C 

DO 58 1*21*41 
11*1-10 

Sill 1 1 l*GED( I l*GEDI I > 

S 1 2 1 1 1 1 *GE C I >*GEC I I 
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S13(II)«FE0(I)*FED(I) 

58 S14IIII*FEIII*FEtIt 

C 

C PERFORM THE INTEGRATION 

C 

STl«SIMPSNCS21*20, HI ) 

ST2»SlMPSNIS22t20,Hl I 
ST3*SIMPSNtS23t2C.Hl) 

T11*PI*SIMPSN(S7»10»H2I 
T12*SIMPSNCSl»20fHl > *S I HP SN < S 8 , 20 , H2 > 
T21*PIPSIMPSN<S9»20»H2) 
T22*SIMPSN!S2t20fHl)*SlMPSN(S10»2C'H2> 
T31«SIMPSNCS3»20»H1 » 

T32«<CEI21l*SEC21l-CEin*SEUM 

T41«PI*SIMPSN(S11»20»H3> 

T42«SIMPSNIS4#20.Hll*SlMPSNI$l2t20»H3l 

T51«PI*SIMPSN(S13»20»H3> 

T52*SIMPSNCS5»20»H11*SIMPSN(S14»20*H3I 

T61«SIMPSN(S6»2C»H1 I 

T62«IFE<41 )*GE(4l)-FEt21t«G£(21)> 

C 

C CALCULATE THE ARBITRARY CONSTANTS. 

C 

A11«U.0D0-0VC2 )/OV(3n*FE(21)*(ST3/ST2» 
A12«C0NS*A(K»*ST1/ST2 

A13*UP*FE(21>»CED(2in/CE(21)-0V(2)9FED(21)/0V(3M 

A21«A11/(A12*A13) 

A1«A21*GE(21 1/SEI21 1 
A1«A1«A1 

B 1 «F e I 21 )*Fe<21 1 
BAl*A21*FE<2l)«GE(21)/S£(21) 

A2«A21*A21 

BA2*A21 

T1«C1*A1*IT11+T12> 

T2«C2*81*IT21+T22 ) 

T3«C3*BA1*T31*T32 
WR I TE ( 6 ♦ 302 I T1,T2»T3 
PC0R*T1*T2-T3 
T4*C4PA2*{T41*T42 > 

T5>C2*1.0DC*(T51*T52 I 
T6«C5*8A2*T61*T62 
WR I TE ( 6 1 302 I T4,T5»T6 
PCLA0«C6*(T4*T5-T6) 

RCOR-PCOR/f PCOR+PCL AO) 

RCLAO-l.ODC-RCOR 

WR I TE ( 6 » 2 1 1 I RC OR t RCL AOt PCOR* PCLAO 
GO TO 70 

81 WRITE(6f212) 

7C CONTINUE 

2 1C F0RMATC1X# 012.51 

211 FORMAT ( 1 IX » 4012. 5) 
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212 FORMAT! 1X» 'NORMAL I ZEO WAVELENGTH « 1* NO POWER CALCULATION' I 

301 FORMAT! lXt 2012.5) 

302 FORMAT! IX, 3012. 5) 

RETURN 
END 

OOUSLE PRECISION FUNCTION S I MPSN! Q »N»H ) 

DOUBLE PRECISION Q!21)»H 

INTEGRATION BY SIMPSON'S RULE 

S I MPSN*Q! l)^4.0D0*Q!2)*Q!N-*l) 

00 1 I *4 » N , 2 

SIMPSN«SIMPSN+2.CD0*Q!1-1 I +4.000*0! I ) 

SIMPSN«SIMPSN*H/3.0D0 
RETURN 
ENO 

SUBROUTINE CHVAL2 ! N » QQ, CHV1 » CHV2 t J ) 

PURPOSE* TO COMPUTE THE CHARACTERISTIC VALUES OF ODD 
ANO EVEN MATHIEU- FUNCTIONS OF POSITIVE OR 
NEGATIVE ' Q* 

INPUT: N-! INTEGER) SPECIFIES THAT CH. VALUES BE 

OBTAINED FOR ORDERS 0 THRU N-l FOR EVEN 
FUNCTIONS ANO FOR ORDERS 1 THRU N-l FOR 
ODD FUNCTIONS 

QQ- ! DOUBLE PRECISION) THE PARAMETER 'O' IN 
MATHIEU'S DIFFERENTIAL EQUATION 
OUPUT* CHV1 - ( D0U3LE PRECISION) AN ARRAY OF LENGTH N 

CONTAINING CH. VALUES OF ODD MATHIEU FUNCTIONS 
OF ORDERS 1 THRU N-l. 

CHVl(N) IS A CUMMY VARIABLE. 

CHV2- ! D0U3LE PRECISION) AN ARRAY OF LENGTH N 
CONTAINING CH. VALUES OF EVEN MATHIEU FUNCTIONS 
OF ORDERS 0 THRU N-l. 

J-l INTEGER) MAXIMUM ORDER UPTO WHICH CH. VALUES 
HAVE BEEN SUCCESSFULLY COMPUTED 
DOUBLE PRECISION CV1 ! 6 , 25 ) »C V2 ! 6, 25 ) ,CHVl !N ) ,CHV2 ! N ) , QQ 
DOUBLE PRECISION QABS * DABS 

131 FORMAT! 'O' »5X, • NOT ALL CH. VALUES AVAILABLE WARNING') 

IFIQQ.LT.O.DO) GO TO 20 

IF!N.GT.l) CALL MFCVAL! N-l »N-1 tQQtCVl »Jl > 

CALL MFCVALIN.N-l ,QQ,CV2» J2) 

IFIJ1.LT.! N-l). OR .J2.LT.N) WRITE!6»1C1 ) 

J*M I NO ! J 1 1 J2-1 ) 

DO 10 I«1»J 

IFIN.GT.l) CHVl(l) «CV1!1, I) 

CHV2! I )*CV2 ! 1 » I ) 

I! CONTINUE 

CHV2! J^l )«CV2! 1 , J^l ) 

RETURN 

20 0A3S*DABS!QQ> 

CALL MFCVAL! N,N-1 »QABS,CV2, J2) 
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IFCN.NE.l) GO TO 25 
CHV2C1 J«CV2C1,1) 

RETURN 

25 CALL HFCVALCN-1 ,N-1 ,QABS,CV1 ,J1 I 

1 F ( Jl.NE.f N-1I.OR.J2.NE.N) WRITE (6, 10H 
J«MIN0(J1,J2-1) 

00 30 1*1 » Jt 2 

CHV2m*CV2fltI I 
CHV1 ( 1 )*CV2( l » I ♦ 1 ) 

CHV2CH-1 )*CV1(1,I) 

I F ( f I +1 ) • LE • J ICHVl ( I +1 ) *C Vl ( 1 • 1 +1 ) 

30 CONTINUE 

IF(M0D{(N-l),2).EQ.0)CHV2(N)*CV2Ci,N) 

RETURN 

ENO 

SUBROUTINE MFCV AL C N , R , QQ« CV » J ) 

C ********** 

INTEGER J*K»KK»L»M,N»R,TYPE 

DOUBLE PRECISION A , C V, DL , OR » DTM » Q, QQ,T » TM»TOL »TOLA 
DOUBLE PRECISION FILLC3) 

DIMENSION C V ( 6 » N ) 

EQUIVALENCE I OL tOR *T ) 

COMMON/MF 1/Q,T0L» TYPE, DUMMY < 4) 

COMMON/MF2/FILL 
T0L*1. 00-13 
IFIN-RI 10,10*20 
10 L * 1 

GO TO 30 
20 L* 2 

3C Q*QQ 

DO 50C K * 1 » N 

J-K 

IFCO) 960,490,40 
40 KK*MIN0(K,4) 

TYPE*Z*M00(L, 2) +MODI K-L*l ,2 1 
C FIRST APPROXIMATION 

GO T0< 100,200*300,400 ),KK 
IOC IFIQ-l. 0001110,140, 140 

112 GO T0( 120,130) ,L 

121 A*1.000-Q-.125D0*Q*Q 

GO TO 420 
13C A«Q*Q 

A*A*I-.50C*.05468750G*A ) 

GO TO 420 

14C IF! Q- 2*0 DO) 150,180,1 80 

15: GO TOC 160,170 ) ,L 

160 A*1 *03300-1 .0 746DC<'0-.G«)89DC*Q j >Q 

GO TO 420 

17: A*.2300-.495D0«'Q-.191D0*Q*Q 

GO TO 420 

1 30 A*-, 2 500-2. ODO*Q^2.:DO*DSQRTC Q) 
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GO TO 420 
200 DL«L 

IFCQ*DL-6.0D0 1 210*350*350 
210 GO T0( 220*230 1 *L 

220 A«4. 01 52 1D0-Q*I.046D0+.C 66785700*0) 

GO TO 420 

230 A« 1.000+1.0500700*0-. 18014300*0*0 

GO TO 420 

300 IFCQ-8.000) 310*350*350 

310 GO T0( 320*330 1 ,L 

320 A«8. 9386700 +.17815600*0-. 0252 13200*0*0 

GO TO 420 

330 A«3. 700 1700 +.953485 00*0-. 047506500*0*0 

GO TO 420 
350 DR«K-1 

A«CVC1 ,K-1 )-DR+4.0DG*DSQRTC0) 

GO TO 420 

400 A*C V C 1 , K- 1 l-C VC 1 » K-2 ) 

A«3.000*A+CVC l»K-3) 

420 IFCQ.GE. 1.000) GO TO 440 
IFCK.NE.il GO TO 430 

425 T0LA«DHAX1C DM1 N1CT0L* DABS 1A1 1,1.00-141 

GO TO 450 

430 T0lA*T0l*DA3S( A! 

GO TO 450 

440 T0LA*T0L*DMAX1 ( QtDABSCA! 1 

445 TOLA*DMAXlC 0M1N1 (TOLA* DABS C A1 * .4D0*DSQRTIQ! 1*1.00-14) 
C CRUDE UPPER AND LOWER BOUNDS 

45C CALL BOUND$CK,A*TOLA»CV,N»M) 

IF(M.NE.O) IFCM-l) 470*910*900 
C ITERATE 

CALL MFITR8(T0LA,CVC1 »K),CVC 2*KI*M) 

1FCM.GT.0) GO TO 920 
C FINAL BOUNDS AND FUNCTIONS* 0 

470 T«CVC l.Kl-TOLA 

CALL TMOFACT,TM,DTH,M) 

IFCM.GT.C) GO TO 940 
CVC3»K)*T 
CVC4,K)«-TM/DTM 
480 T*CVI l »K l+TOLA 

CALL TMOFA(T,TM,OTM,M) 

IFCM.GT.C) GO TO 950 
CVC5*K)«T 
CVC6,K)*-TM/0TM 
GO TO 5C0 

C 0 EQUALS ZERO 

49C CVC1 ,K)«CK-L+1 1**2 

C VC 2 »K 1*0. 000 
CVC 3*K )»CV ( 1 *K 1 
CVC4,K)«0.0D0 
CVC5,K)«CVC1,K) 
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CV!6,K)«0.000 
500 CONTINUE 
550 RETURN 

C PRINT ERROR MESSAGES 

90? WRITE! 6»901 1 K 

9C1 FORMAT! *0* * 'CRUDE BOUNDS CANNOT'** BE LOCATEO* NO OUTPUT', 

C ' FOR K« • » I 2 ) 

GO TO 930 

910 WRITE!6,911) K 

911 FORMAT! »D' , 'ERROR IN SUBPROGRAM TMOFA* VIA SUBPROGRAM 

C BOUNDS, NO OUTPUT, FOR K»'*I2I 

GO TO 930 

92C WRITE! 6,9211 K 

921 FORMAT! 'O' , 'ERROR IN SUBPROGRAM, TMOFA, VIA SUBPROGRAM, 

C MFITR8, NO OUTPUT, FOR K«*,I2I 

930 J«J-l 

GO TO 550 

9*0 WR I TE!6»9*1) K 

9*1 FORMAT! 'O' , 'ERROR IN SUBPROGRAM, TMOFA, NO LOWER BOUND, 

C FOR K*' » 1 2 I 

C V 1 3 , K )*0 • ODO 
CV!*,KI«0.00 
GO TO *30 

950 WRITE(6,95ll K 

951 FORMAT! 'O' , 'ERROR IN SUBPROGRAM, TMOFA, NO UPPER BOUND, 

C FOR K«' ,12 ) 

CVI5»K) = D*00 
CV! 6,K 1*0. DO 
GO TO 5C0 

96? WRITE(6,961 ) 

961 F ORMAT ( 2CH3Q GIVEN NEGAT I VELY, , 20H USED ABSOLUTE VALUE) 

Q*-Q 

GO TO *0 
END 

SUBROUTINE BOUNDS ( K , APPROX , TOL A ,CV ,N, MM I 
INTEGER K,KA,M,MM,N 

DOUBLE PRECISION A , APPROX , AO , A1 , CV , DTH , DO , D1 , 0, TM ,TOL A 
DIMENSION CV ! 6 , N ) 

C0MM0N/MF1/Q, DUMMY (7> 

C OMMON/MF 2 / AO » A * A 1 
KA«0 

IF(K.EQ.l) GO TO 2C 
IF! APPROX-CV! 1,K-1 ) ) 1C, 10, 20 
ID AO*CV( 1 ,K-1 I^l.ODO 

GO TO 3C 
20 AO-APPROX 

3D CALL TMOFA! AO, TM,DTM,M) 

IFIM.GT.OI GO TO 250 
DO--TM/DTM 
I F ! 00 ) 100,300,50 
C AO IS LOWER BOUND, 



c SEARCH FOR UPPER SOUND 

50 A1 *AO*DO+« 100 

CALL TM0FAIA1,TM,DTM,M> 
l F ( M. GT.O ) GO TO 250 
Dl — TM/DTM 
IFCD1I 200 » 350 * 60 
60 A0*A1 

D0*01 
X A«KA*1 

I F CXA-41 50*400*400 

C Al IS UPPER BOUND* SEARCH FOR LOWER SOUND 

ICO A 1 *A0 

01*00 

A0*DMAXHA1 + D1-.1D0*-2.0DC*Q> 

IF(K.EQ.l) GO TO 110 
IFUO-CVI l.X-1 1> 150*150*110 
110 CALL TMOFA(AO»TM*DTM»M) 

IFIM.GT.OI GO TO 250 
DC*-TM/DTM 
I F ( DO ) 120 »30C f 200 
120 X A*K A+l 

I F ( X A-4 I 100*400*400 
150 XA«XA+1 

IFIXA-41 160*400*400 
160 AC*Al4-0MAXl( TOLA* DABS CD1 ) ) 

GO TO 30 

2C: A*.5DC*C A0+D0*AltDl I 

IFIA.LE.AO.OR.A.GE.A1) A* . 5 DC* ( A0-* A 1 I 
250 MM*M 

RETURN 

30 0 CVIl»XI*AO 
310 CV(2»X)*0.00 

M*-l 

GO TO 250 
35C CV(1*X)*A1 

GO TO 310 
400 M*2 

GO TO 250 
END 

SUBROUTINE MFITR8(T0LA,CV*DCV»MMI 
INTEGER M » MM »N 

OOUSLE PRECISION A * AO * Al * A2 * CV » D» DC V*DTM * TM »TOL A 
LOGICAL LAST 
COMMON/MF2/AO»A * Al 
N*0 

LAST*. FALSE. 

50 N*N«-1 

CALL TMOFAI A,TM,DTM*M) 

IFIM.GT.OI GO TO 4CC 
D*-TM/OTM 

C IS TOLERANCE MET 
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IFIN.EQ.AO.OR.A-AC.LE.TOLA.OR.Al-A.LE.TOLA.OR. 

C DA8SC DI.LT.TOLA)LAST«.TRUE. 

I F f 0 1 110.109*120 
10C CV*A 

ocv«o.oo 

GO TO 320 

C REPLACE UPPER BOUND BY A 

110 A 1 «A 

GO TO 200 

C REPLACE LOWER BOUNO BY A 

1 20 AO«A 
2 DO A2«A*0 

I FILASTl GO TO 300 

IF(A2. GT.A0.AN0.A2.LT.An GO TO 250 
A«.5D0*( A0+A1I 
GO TO 50 
250 A*A2 

GO TO 50 

300 IF( A2.LE.A0.0R.A2.GE.A1) GO TO 35C 
CALL TM0FA(A2,TM*DTM,M> 

IF(M.GT.O) GO TO 40C 
0«-TH/DTM 
CV«A2 
310 DC V*0 

32C MM*M 

RETURN 
350 C V* A 

GO TO 310 
400 C V = 9 • DO 

DCV*0 • DO 
GO TO 320 
END 

SUBROUTINE TMOF A ( AL FA ,TM , OTM * NO » 

INTEGER K.KK,KT,L*MF,MC»M1 »M2S«ND»TYPE 
OOUBLE PRECISION A , AA , ALF A * B ,DG ,DTH, DTVP E « 

C F»FL»G*HC200l*HPtQ»QlNV» 

C 01 »02 »T »TM»TOL*TT »V 

COMMON G( 200.2) ,DGt 2C0.2 I ♦ A A , A ( 3 1 • B( 3 ) t DTYP E t OINV » 
C 01,Q2»T,TT,K,L»KK»KT 

C0MM0N/MF1/Q,T0L.TYP5»M1 ,M0»M2Sf MF 
EQUIVALENCE C HID » GC 1 * U 1 1 1 01 1 HP I . C 02 » F 1 
DATA FL/1. 00*30/ 

C STATEMENT FUNCTION 

V(KI«C AA-D6LEI FLOAT! Kl)«2)/Q 
ND»0 
KT«0 
AA«ALFA 
DTYPE*TVPE 
0INV*1 .ODC/O 
DO 10 L*1 *2 
DO 5 K»1»2C: 
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GCX»Ll»O.DO 
DGCK»LI*0.DC 
5 CONTINUE 

1C CONTINUE 

I F C MODI TYPE 1 2 I ) 20t30#20 
20 M0»3 

GO TO 40 
30 M0«TYPE+2 

40 K*.5D0+DSQRTC DMAX1 {3.CDO*Q + AA»0.00 I » 

H2S«MINO(2*K + MO+4 t 398*MOO(MO >2 1 1 

C EVALUATION OF THE TAIL OF A CONTINUED FRACTION 

All )«1.0D0 
A ( 2 l*VI H2S+2 I 
Btl )*VCH2S) 

BC2I«AC2I*3(1 I - 1 .000 
Ql«A(2)/B(2) 

DO 50 K*1 » 200 
HF-M2S+2+2** 

T«V(NF» 

A C 3 I *T*A C 2 I- A < 1 ) 

BC3I*T*BC2»-BI1I 
Q2«A(3)/B< 31 

IFI0ABSCQ1-Q2I.lt. TOD GO TO 70 
Q1*Q2 
A C 1 I ■ A C 2 I 
A ( 2 » « A C 3 » 

B « 1 > « B C 2 > 

BC2I*3I3) 

50 CONTINUE 

KT«1 

70 T*1 .0D0/T 

TT»-T«T*Q1NV 
L*HF-M2S 
00 80 K*2 ♦ L # 2 
T*1.0D0/CVCMF-KI-T I 
TT«T*T*CTT-QINV I 
80 CONTINUE 

KK*M2S/2*1 
IF(KT.EO.l) Q2«T 
G(KK»2I*.5DC*(Q2*T> 

DGC KK.2 )*TT 
C STAGE 1 

GC2» 1)«1.0D0 
DO 140 K«H0»M2S*2 
KK*K / 2*1 

IFCK.LT.5) 1FCK-3) 1CC»11C»120 
G ( KK » 1 ) *V dC-21-l.CDC/GCKK-ltl ) 

DGCKKtl |«Q1NV + DGCKK-Itl )/G(KK-l»l 1**2 
GO TO 130 
ie: G(2«n«vco) 

DGC2t ll-QINV 



CO TO 130 

110 G(2,l »*V( 1H-0TYPE-2.0D0 

06(2*1 »«QINV 
GO TO 130 

12C G(3»l l«VI 21 ♦COTVPE-2 •00I/G12 *1 > 

0G(3»1 )«QlNV+( 2.00-DTYPE l*DG(2 *1 )/G(2* 1 1**2 
IF(TYPE.EQ.2» 6(2*11*0.00 
13c IF(DA5S(G(KK, in.LT.l.DOl GO TO 200 
14C CONTINUE 
C BACKTRACK 

TH«G(KK,2 l-G(KK *1 I 
DTH*DG(KK,2I-DG(KK,1 ) 

H1«H2S 

KT«H2S-M0 

DO 180 L* 2 * KT» 2 

K*M2S-L 

KK*K/2+l 

G(KK,2»«1.00/(V(K»-G(KK*l*2) ) 
0G(KK»2l*-G(KK»2)**2*(QINV“DG(KK'»l,2ll 
I F ( K-2 I 150*150*160 
150 G(2,2>*2.0D0*G(2,2 ) 

DG( 2 * 2 )*2. DO*DG (2*2) 

160 T*G(KK,2»-G(KK,1I 

IF(DABS(T)-0A3S(TMI I 170»18:,18C 
170 TM*T 

OTM*DG( KK, 2 |-DG(KK * 1 ) 

H 1 *K 

180 CONTINUE 
GO TO 320 
C STAGE 2 

200 M 1 *K 

K*M2S 
KK*K/2+l 

210 I F(K.EQ.Hl ) IF(K-2» 3CC»300.31C 

K*K-2 
KK*KK-1 

T*V(KI-G(KK+1,2 I 
!F(OA8S(TI-1.00 » 250,220*220 
22- G(KK,2I*1.000/T 

DG(KK,2I*(DG(KK+1»2)-QINV>/T**2 
GO TO 21C 
C STAGE 3 

250 IF(K.EQ.Ml) IF(TI 220*290,220 

HP*DG(KK*1 ,2 I-QINV 
260 G(KK*2 I *FL 

H(KK)*T 
K*K-2 
KK*KK-1 
F* V( K I *T- 1 • DC 

IF(K.EQ.H1 ) I F ( F » 280,290,280 
IF(0A3S(FI-0A3S(T) I 270*280,280 
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27D 

2 80 

2 90 
C 

300 

31C 

320 

C 

C 

C 


101 

102 

C 

C 

C 


5 

C 

C 

c 


HP»HP/T**2-QINV 


T«F/T 
GO TO 260 
G!KK»2 >*T/F 

DGtKKt2l«<HP-QINV6T*T>/F**2 
GO TO 210 


NO* 1 

GO TO 320 

CHAINING M EQUALS 2 
G! 2»2I*2.00*G!2 #2 I 


DG!2#2l*2.D0*DG!2»2 ) 

TH*GlKK»2)-G(KKtl > 

DTM*DG!KK#2)-DG!KK»1 ) 

RETURN 

END 

SUBROUTINE EXP AND! QD» PNC » R #CV»NORM »CD»N I 

PURPOSE: TO GET EXPANDING COEFFICIENTS FROM ROUTINE 

COEF. TERMINATE THE TERMS FOR REQD. ACCUTACY 
ANO DO THE NORMALIZATION 


DIMENSION CD!25> „ 

DOUBLE PRECISION A » CV » QD » Q» TOL »T « ABt ERR »DA8 S » 

C CD»SUH»Tl * D SORT * SUMl 

INTEGER R#FNC» TYPE# CASE# NORM 
COMMON DUM1 1 1600 I » A tT »DUM2 1 6 I » A8I 200 l 
COMMON /MF1/Q»T0L#TYPE»M1»MD#M25»MF 
FORMAT! *0* » *THE * OF EXPANDING COEFFICIENT REQD. 

C THAN 25* »5Xt* WARNING* I 
FORMAT! *0* . 'ERROR IN SUBPROGRAM* » ' TMOF A VIA COEF. 

C ARGUMENTS NO OUTPUT* I 
T0L*l.0D-13 

TO TEST THE LAST VALUE OF ARRAY CD ERR IS USED 


IS MORE 
VERIFY 


ERR*1.0D-20 


Q = QD 

TYPE*2*MOO!FNCf2)+MOD!R*2 I 

FOR NEGATIVE 0 AND ODD ORDERSt EXP. COEFFS. FOR EV:N AND 
ODD FUNCTIONS ARE INTERCHANGED. 

IF!Q.LT.O.ODO.AND.MOD!R»2I.EQ.l ITYPE*2*M0D1 !FNC-1)#2I*H0D! 
M*0 


A*C V 

Q*DABS!Q) 

CALL COEF! M I 
I F I M.EQ.O ) GO TO 5 
MR 1 TE ! 6 » 102 ) 

RETURN 

TYPE*2*M0D! FNC# 2 I ♦MOD! R 1 2 I 

THE^COEFFICIENTS PASSED THRU COMMON ARRAY A8 IN DOUBLE 
PRECISION IS GIVEN TO AN ARRAY CD OF LENGTH 25 FOR 
FURTHER PROCESS ION 
DO 10 1*1 #25 
CD! I I *AB ! I > 


69 


IF(DA8S(C0m).LT,l.D-30> CO <11*0.00 
10 CONTINUE 

IF<CD(25).GT.ERR) WRITE<6»101) 

Nl-R/2+1 
00 20 I *N1 f 25 

IF(CO( D.EQ.O.OO) GO TO 25 
N* I 

20 CONTINUE 

C NORMALISING THE CODES. PRESENTLY IN NEUTRAL NORM 

25 SUM*0.000 

I F( NORM. EO. 1 I GO TO 14C 
C GETTING STRATTON NORMALISATION FACTOR 

IF(QO.LT.O.OO) GO TO 91 
GO TO <40i40»60»8C)*CASE 
40 DO 50 J* 1 1 N 

SUM«SUM^CO< J> 

50 CONTINUE 

GO TO 100 
60 00 70 J*ltN 

SUM«SUM+CD< J»*0BLE< FLOAT! 2*1 J-i 1 1 I 
70 CONTINUE 

GO TO 100 
80 DO 90 J* 1 « N 

SUM*SUM»CD(JI$D3LE (FLOAT < 2® J-l ) I 

90 CONTINUE 
GO TO ICO 

C GOT NEGATIVE 0 STRATTON NORMALISATION FACTOR IS DIFFERENT 

91 T 1 *- 1 . ODO 
IF(M0D<R/2»2).E0.1 )T1=-T1 
GO T0< 92t92,94 t 96) ,CASE 

92 DO 93 J* 1 » N 
T1«-T1 

SUM*SUM*T1«CD< J ) 

93 CONTINUE 
GO TO 100 

94 DO 95 J* 1 1 N 
T 1 *-T 1 

SUM*SUM*CD< J)*T1*03LE(FL0ATC2*| J-l I 1 1 

95 CONTINUE 
GO TO 100 

96 00 97 J* 1 « N 
T 1 *-Tl 

SUM = SUM + CD(J)FT1*08LE< FLOAT < 2*J-1I I 

97 CONTINUE 

ICO I F < NORM. EO. 2 I GO TO 120 

C GETTING I NCE ' S NORMALISATION FACTOR 

SUM1*0. DO 
DO 110 J*1 1 N 
SUM1*SUM1«-C0< J >®CD< J) 
li: CONTINUE 

IF(FNC.EQ.2.AND.M00<R»2).E0.C') SUM 1 *SUM1 +C0 < 1 ) *C0 ( 1 ) 



non*-* 
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C 

12C 
30 


140 

is: 


c 


6C 


7C 


80 

90 


no 

no 


DSIGNCSUMltSUMl 
BY NORMALISING FACTOR FOR 2 


TYPE I C ASE*3 I 
COEF RETURNS 


SUM1-0SQRTCSUM1I 
IF(N0RM.EQ.3I SUM 
01 VI OE ALL COEFS. 

DO 130 1*1 tN 
cdui*co< ii/sum 
CONTINUE 

FOR MATHIEU FUNCTIONS OF SE2N+2 
BE B2» 64 ETC* BUT THE ROUTINE 
THIS IS TO BE DROPPED* 

IFICASE.NE.3I RETURN 

DO 150 1*2 »N 

C0C 1-1 l*CD( I I 

CONTINUE 

CDINl*0.D0 

RETURN 

ENO 

DOUBLE PRECISION A,AB*FL*Gt M C20^l*Q*T*TOL.V 
COMMON GC 20C *2 It DUM1 C 800 » * A *T#K»KA* KB»KR *MM 
COMMON /MFl/QtTOLtTYPEtMl ,MCtM2StMF 
EQUIVALENCE IH< 1 1 tGC 1 tl 1 1 
DATA FLt V2/1.0+30* 1.0-15/ 

STATEMENT FUNCTION 
VCK)*(A-DBLE(FLOATIK) I **2 I / Q 
CALL TMOFA(AiT,T,M) 

I F < M.NE.O I GO TO 30C 
DO 60 K* 1 t 200 
ABI 10*0.00 
CONTINUE 
KA = Hl-M0«-2 
DO 90 K*2 t KAt 2 
KK*CMl-Kl/2«-l 
IFCK-21 70t70t9C 
abikki*i.do 
GO TO 90 

A3<KK)*AB<KK«-l l/G(KK + l tl I 

CONTINUE 

KA*0 

00 130 K*M1 . M2S t 2 
KK*K/2^1 


HL s K 

IFIGIKKt2l.EQ.FL) GO TO ICC 

ABIKKI*ABCKK-1)4G(KK»2) 

GO TO 11C 

T* AB ( KK-2 ) _ _ 

1FCK.EQ.4.AN0.M1 .EQ.2 I T*T+T 
AB(KKI*T/IV(K-2)*HIKKI-1.D0) 

IFC DABS! AB(KK)) .GE.l . 0-17 I KA*C 

IF(KA.EQ.5I GO TO 260 
KA*KA*1 


£ 3 ONLY 


COEFS. SHOULD 
A 60*0 ALSO. 


tV2 

,MLtAB(2C0l 
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130 CONTINUE 

T*DLOG(DABS(AS(KK) » /V2 I /OLOG 1 1 .DO/DABS I GC KK 1 2 > I ) 

KA«2*I DINT IT I 
ML«KA*2«-M2$ 

I FI ML* GT • 399 1 GO TO 403 

KB*KA*2+MF 

T-1.00/VIK8I 

KK*HF-M2S 

00 150 K«2»KK,2 

T«1*03/IVIK9-K)-T> 

150 CONTINUE 
KK«ML/2+l 
GIKK,2 )*T 
00 203 K*2«KA»2 
KK* I ML-K I /2+ l 

G < KK » 2 )*1*D0/(V(ML-K)-GIKK*1 ,2 I ) 

23C CONTINUE 
KA-M2S+2 

DO 253 K*K A , ML » 2 

KK»K/2+l 

A B I KK ) * AS I KK- 1 l«GIKK,2 > 

250 CONTINUE 
C NEUTRAL NORMALIZATION 

263 T«ABI1) 

MM*MODI TYPE * 2 } 

KA*MM«-2 

DO 283 K*K A , ML » 2 
KK*K/2+l 

IFIDABSITI-OABSIABIKK) ) ) 270,280,280 
2 70 T s A3 1 KK ) 

MM*K 

280 CONTINUE 

DO 29C K* 1 » KK 
AS I K ) *AB I K ) / T 
293 CONTINUE 
33C RETURN 
430 M«-l 

GO TO 300 
END 

00U3LE PRECISION FUNCTION ANGMFC I QD * FNC ,R , XO , DSA l V , AR , N I 
C PURPOSE: TO COMPUTE A PERIODIC MATHIEU FUNCTION, 

C ODD OR EVEN TYPE OR ITS DERIVATIVE 

DOUBLE PRECISION PC , PS , DP C , DPS 
EXTERNAL PC , PS , DPC , DPS 
DIMENSION ARI25),ABI25) 

INTEGER FNC, R, DERI V, TYPE, CASEtP 

OOU3LE PRECISION AR , A 3 , XD , X , T 1 , SUM , DCOS »0S I N , QD 

COMMON/NTERM/NL IMIT 

COMMON/ANG/AB , X ,P 

NLIMITsN 

QS*QD 



17 


1 

C 

c 

20 


30 

35 


40 


5 0 


60 

70 


8C 


9C 

10C 

110 

120 

13C 

1*J 


C 


X«XD 

TVPE*2RMODCFNC,2 I ♦MODI R * 2 ) 

CASE-TYPE+i 

00 1 I«1»N 

ABC II«ARC II 

CONTINUE 

FOR NEGATIVE Q IN ALL SUMMATIONS 
A MINUS SIGN 
I F C QS I 20 * R0» 35 
Tl — l.ODO 

IFCCASE.EQ.3IT1-1.00C 

IFCM0DCR/2.2I.EQ.I IT1«-T1 

00 30 I«1»N 
Tl— T1 

ABC 1 I*T1*AB( 1 I 

CONTINUE 

P»-l 

IFCCASE.E0.1> P*"2 
IFCCASE.EQ.3I P«0 

1 F COER IV.EQ.l I GO TO 60 
GO T0C40»4C*50*SCI .CASE 
CALL S I GMAC PC * SUM I 
ANGMFC-SUM 

RETURN 

CALL S I GMAC PS » SUM I 

ANGMFC*SUM 

RETURN 

GO TOC TO *70 * 8C * 80 ) .CASE 
CALL S I GMAC OPC * SUM ) 

ANGMFC«SUM 

RETURN 

CALL S I GMA C OPS » SUM ) 

ANGMFC-SUM 

RETURN 

IF10ERIV.EQ.DG0 TO 120 
GO TOCIOO, 100*110, HClfCASE 
ANGMFC*0C0S COBLECFLOATCR) !*X) 


ALTERNATE TERMS HAVE 


RETURN 

angmfc*dsincdblecfloatcri »* x I 
RETURN , _ 

GO TOC 130*130, 140*140I»CAS£ 
ANGMFC — R«OSINC DSL EC FLOAT C R*> l*X 


I 


RETURN 

angmfc*r*dcos 


COBLE C FLOAT C R I l*X I 


RETURN 

END 

DOUBLE PRECISION FUNCTION PCCKI 

INTEGER p*x 

DOUBLE PRECISION AB l 2 5 ) ♦ X * DCOS 


COMMON/ AN G/ AB » X »P 

EVALUATES ONE TERM OF THs 


EVEN PERIODIC SOLUTION 


no on 
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P C« AB( K)*DCOS ! 03L EC FLOAT !2»K*P) )*X) 

RETURN 

END 

OOUBLE PRECISION FUNCTION PS!K) 

INTEGER PtK 

OOUBLE PRECISION AB! 25 ) *X ,DS IN 
COMMON/ ANG/AB»X fP 

C EVALUATES ONE TERM OF THR ODO PERIODIC SOLUTION 

PS*AB!K)«DSIN!DBLE!FL0ATI2*K+P ) )*X) 

RETURN 
ENO 

DOUBLE PRECISION FUNCTION DPCIRI 
INTEGER P'K 

OOUBLE PRECISION AS I 25 ) , X .T , DS IN 
COMMON/ANG/ABf XtP 

EVALUATES ONE TERM OF THE DERIVATIVE OF THE EVEN PERIODIC 
MATHIEU FUNCTION. 

T*2*K+P 

DPC*-AB< K)*T*DS IN! T«X) 

RETURN 
END 

DOUBLE PRECISION FUNCTION DPSIKI 
INTEGER P,K 

DOUBLE PRECISION ABC 25) »X »T »DCOS 
COMMON/ANG/AB»XtP 

EVALUATES ONE TERM OF THE DERIVATIVE OF THE ODD PERIOOIC 
MATHIEU FUNCTION. 

T«2*K*P 

DPS*AB!K)*T*DCOS!T*X ) 

RETURN 
END 

SUBROUTINE S I GM A( DUM , SUM ) 

C PURPOSE: TO SUM N TERMS CSPECIFIED BY THE COMMON 

C COMMON BLOCK N TERM) OF A FUNCTION. 

OOUBLE PRECISION DUM, SUM, ERR ,T1 .TERM, DABS 
COMMON/NTERM/NL IM IT 

1C1 FORMAT! 'C**' CONVERGENCE NOT TO S AT I SF ACT I ON * 1 1 0 X » • WAR N I NG * ) 

ERR«1.0D-13 
T1*DUM! I ) 

SUMsTl 
N*NL I M I T 

IF! NLIMIT.GE.22 )N*22 
I M I N* 5 
DO 10 I«2 »N 
TERM*DUM ( I ) 

SUM=SUM*TERM 

TERM-DABSCTERM) 

IF! I.LT.IMIN) GO TO 1C 

IF10ABS! SUM).GE.ERR*TERM) RETURN 
I 0 CONTINUE 

20 IF! I.EQ.22 )WRITE!6»10l ) 
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C 

13 


20 


C 

c 

30 


40 


C 

53 


60 


RETURN 
END 
SUBROUTINE 


FACTOR* FNC,R *00* AB»N, PS I 


PURPOSE* 


TO COMPUTE THE NORMALIZATION FACTOR»P2N»P2N« , l • 
S2N+1*0R S2N42 FOR MODIFIED MATHIEU FUNCTIONS 
OF POSITIVE *Q' AND P2N* *P2N*1 • *S2N*1 • t 
OR S2 N+2 * FOR THOSE OF NEGATIVE *0* 

DIMENSION AS* 25 ) 

INTEGER FNCfRtCASE 

DOUBLE PRECISION AB , PS, PSP , DABS .OSQRT *RQ* ODEV t 
c SUM 1 ,SUM 2 »T 1 ,T 2 *QD 

RQ*DSORTC DABS CQDI I 
CASE« 2 *M 0 DCFNC» 2 I 4 M 0 D*R. 2)+1 
IF* QD*LT*O.ODO )CASE*CASE 4 4 .' 

ODEV*1*ODO 

IF(M0D*R/2»2I«cQ.1 I0DEV*-1*0D0 

SUM1*0.0D0 

SUM2«0.CD0 

T1--1 .ODO 

GO TO < 10*30, 50*70* 10, 70*50, 30I»CASc 

FOR ALL Q ANO EVEN ORDER P2N AND P2N* FNC*2 

00 20 1*1, N 
SUMI*SUM1*AB( I ) 


T1*-T1 

SUM2*SUM2^T1«ABC I I 
CONTINUE 

PS*SUM1*SUM2/A9< 1 1 

PSP*P$*ODEV 

IFCQD.LT.O.ODOIPS*PSP 

FOR POSITIVE Q AND ODD ORDERS P2N+1, P2N+1* IF FNC*2 
NEG Q AND ODD ORDER IF FNC“1 
DO 40 I * 1 , N 
Tl—Ti 

SUH1*SUM1+ABIII 

SUM2«SUH2«-T1*A9* I l*DBLE * FLOAT* 2*1-1 1 > 


CONTINUE 

PS*SUMl*SUM2/*RQFABtll) 
PSP*PS*OOEV 
1 F ( 0D« LT .0 • DO I PS-PSP 
RETURN 

FOR ALL 0 AND EVEN ORDER 


IF FNC* 1 


S2N + 2 » S2N*2* 


Tl*l .ODO 
DO 60 I * 1 ♦ N 

T2*AB( I »*DBLE* FLOAT *2*1 1 1 

SUM1-SUM1+T2 

Tl— T1 

SUM2*SUM2*T1*T2 

CONTINUE 

PS«SUM1*SUM2/(RQ*RQ*AB(1I 1 

PSP«PS*ODEV 
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IF<Q0.LT.0.CD0IPS*PSP 

RETURN 

C FOR POSITIVE Q AND OOO ORDER IF FNC«1 

C FOR NE6ATIVE Q ANO ODD ORDERS S2N*l» S2N*1* IF FNC-2 

70 00 80 I«1.N 

SUH1*SUMI*ABI l l*DBLE( FLOAT (2*1-1 » I 
Tl — T1 

SUH2*SUM2*Ti*ABU I 
eo CONTINUE 

PS*$UH1*SUM2/(R Q*A 6(1 ) > 

PSP«PS*ODEV 

I FI QD. LT ,0. 000 I PS «PSP 

RETURN 

END 

SUBROUTINE STOR E( QD» X I » NMAX > 

C PURPOSE: TO COMPUTE AND STORE THE VALUES OF BESSEL 

C FUNCTIONS AND DERIVATIVES REQUIRED IN HATHIEU 

C FUNCTION CALCULATION. 

00U3LE PRECISION QD » QABS ♦ RQ, DABS* DSQRT, DEXP ,BS I VI , BS l V2 , 

• BS2V2»DB$lVl,DBSlV2,DBS2V2tXl ,Vl»V2 
COMMON/LOCAL/DUMMY1 (41 ,V1 » V2 * DUMHY2 ( SD ) 

COMHON/R ADI AL/SS1 VI ( 25 1 »3S1 V2 ( 25 ) » 3S2V21 25 ) * DBS1 VI ( 25 ) t 

* DBSIV2I25I,0BS2V2(25I 
N*NHAX*3 

IFIN.GE.25 )N«25 
Nl-N-1 

QABS*DA3SIQD) 

RQ*DSQRT ( QABS } 

Vl*RQ*DtXP<-XI » 

V2*QABS/V1 

IFC QD.LT.O.ODC » GO TO 20 
CALL BESSELfltVltBSIVltNI 
CALL BESSEL(ltV2»BSlV2»N) 

CALL BESSEL(2«V2tBS2V2»N) 

DBSIV1 (I) — BS1VK2 I 
DBS 1 V2 ( 1)«-8S1V2(2 I 
0BS2V2C 1 )»-BS2V2(2) 

00 10 I«2»N1 

DBS1V1 ( I )«( BS1VK 1-1 )-BSlVH 1*1 1 1*0.500 
DBS 1V2 ( I ) * ( BS1 V2 ( I - 1 I-BS1V2I 1*11 1*0.500 
D3S2V2 ( I ) * ( BS 2V2 ( 1-1 1-3S2V2! 1*1 1 )*G.500 
10 CONTINUE 

RETURN 

2C CALL 3SL2(l«VltBSlVl*N) 

CALL 3SL2(ltV2t8SlV2tN) 

CALL 8SL2I2.V2»8S2V2»N> 

DBS1V1C1 I-BS1V1 (21 
DBS1V2(1)*BS1V2(2I 
DBS2V2C1I—8S2V2I 2 I 
DO 30 I *2 1 N1 

DBS1V1 (I)«( BSIVlf 1-1 I+BS1VH 1*1 1 1*0.500 
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DBS1V2C I WBS1V2I I-1I+3SIV2C l*\ » 1*0. 5DC 
03S2V2 ( II*“CBS2V2C 1-1I+6S2V2C1*1I 1*0*5 DO 
39 CONTINUE 
RETURN 
ENO 

SUBROUTINE 8ESSELC SOL.U*BSJY»N> 

INTEGER NtNN.SOL 

DOUBLE PRECISION BSJYCN),U 

NN*N-I _ „„ 

I F ( U* EQ« 0* 00* AND# SOL • EQ» 2 ) GO TO 80 

I F C U*GE • 8» 09 I GO TO 30 
GO TOC 10*201 »SOL 
CALL JCJICUiBSJY) 

GO TO 40 

CALL Y0Y1 ( U*B$ JY) 

GO TO 49 

CALL LUKECUfSOL.BSJY* 

IFCN.LT.2 > GO TO 100 
GO TOC 50 t60 I »SOL 
CALL JNSCBSJY,U,N) 

GO TO 100 

RECURRENCE FORMULAR 

g5 jY^ k+ 1 I *2 •DC*D3LE C FLO AT C K- II l*BS JY C K I / U-BS J V C K* 1 > 
CONTINUE 
GO TO 100 
NN*NN+1 
00 90 K*1 1 NN 
BSJY(KI«-1 .9+37 
CONTINUE 
RETURN 
END 

SUBROUTINE JOJICXtBJ) 

DOUBLE PRECISION BJC2l#TC51tX 
T C 1 1 * X/ 2 • DO 
BJC1)«1*00 
BJC 2»«TC 1 » 

T C 2 l*-T C 1 1**2 
T C 3 I “1 • 00 
TC4»«1.D0 

1C TC4)«TC4)*TC2I/T(3)**2 

BJCII«BJCIH > TC4! 

TC5l«TC4l*TCll/CTC3l*l.D9l 

JFCDMAX1I0ABSCTI4) » « DABS C T C 5 H I .L T . 1 . 0-1 5 I RETURN 
TC3I*TC3I*1.D0 

GO TO 10 
ENO 

SUBROUTINE YCY1 (XtBYI 
DOUBLE PRECISION T 1 1 C I » X , BY t 2 I 
T( 1 )*X/2.D0 


10 

20 

30 

40 

50 

C 

60 

70 

30 


90 

ic: 


o o 
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TI 2 l«-Tin «2 
BY! 11*1.00 
SYC 2 l*T 1 1 1 
T C 7 1 *0 • 00 
T< 10 l*-T 1 1 ) 

T ( 31*0.00 
T 1 4 1 *0 • 00 
TI5I-1.00 

1C T ( 3 >*T C 3 1 ♦ 1 • 0C 

TI4I*T(4I+1.D0/T(3) 

T(5I*TI5I*T(2)/TC3I**2 
BY(ll*BVtl)*TI51 
TI6I--TI 5)*T(4> 

T( 7 l*Tf 7 >*T( 6 ) 

T I 8 1 *T I 5 I *T ( 1 1 / f T I 3 1*1 • DC I 
BYt2l*9Y(2)*T(B> 

T<9)«-Tt 8l*t2«DC*TUI*l.DG/(TI3>*l.DCI i 
TI1C»«TU0I*TI*I 

IF( DM AX1(DABS(T(6)),DABS(T(9H).GE. 1.0-15) GO TO 10 
TI2J*. 577215664 90 153286 00 ♦DLOGCTI 1 1 1 
BY! 1 I *.6366 197723675 8 13400* I BYI1I *T C 2 » «-T I 7 > ) 

BYI2 )*. 63 66 1977 2 367 5 31 3400* I BY(2)*T< 2 )- 1 . DC /X ) *T ( 1 0 1 / 

C 3.141592653589793200 
RETURN 
END 

SUBROUTINE LUKE ( U , K I ND» SS JY I 
INTEGER K • K I NO 

DOUBLE PRECISION A 1 1 9 » , 9 { 1 9 ) , CS , C ( 1 9 ) , 01 1 9 > ,G ( 3 ) , BS J Y ( 2 ) , 
C R(2),S(2),SN,T,U,X 

WARNING - THE FOLLOWING DATA STATEMENTS ARE NOT IN ASA 
STANDARD FORTRAN 
DATA A/. 9995950647696728741600. 

* -.538079561396069130-3. 

* -.131796771233615700-3. 

* .151422497C48644D-5. 

* .158468617920630-6. 

* -.856C695539460-8. 

* -.295723433550-9, 

* .65735562540-10, 

* -.223749703D-11 , 

* - . 4482 1140 D- 12, 

* .69548270-13, 

* -.1513400-14, 

* -.924220-15, 

* .155580-15, 

* -.4760-17, 

* -.2740-17, 

* .610-18, 

* -.40-19, 

* -.10-19/ 

DATA 3/-. 7769355694205321360-2, 
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* 

-*7748032309654476700-2 » 

* 

•25365411654307960-4, 

* 

•394273598399711 D- 5, 

* 

-.107234982991290-6, 

* 

-.7213897993280-8, 

* 

•737646028930-9, 

* 

•150687811D-11* 

* 

-.5745895370-11, 

* 

•459965740-12, 

* 

.22703230-13, 

a 

-•8878900-14, 

* 

•744970-15, 

* 

.58470-16, 

* 

-•24100-16, 

* 

.2650-17, 

* 

.130-18, 

* 

-•100-18, 

* 

.20-19/ 

DATA 

C/ 1.000 67753 5 8659 134623 400, 

* 

.901007251959081830-3, 

* 

.221724349185994540-3 , 

* 

-.1965759463191040-5, 

* 

-.20889531143270-6, 


.10281443508940-7, 

* 

.37597C547390-9, 

* 

-.7638891358D-10, 

« 

. 2387346700-11 , 

« 

.518254890-12, 

» 

-.76939690-13, 

* 

.1440080-14, 

* 

.1032940-14, 

« 

-.168210-15, 

4 

.4590-17, 

4 

.3020-17, 

4 

-.650-18, 

4 

.40-19, 

4 

.10-19/ 

DATA 

0/. 2 3 37632 993 6 2 85 8 03 2 8 0-1 , 

4 

.23346801223545575330-1, 

4 

-.3576010 5909013820-4, 

4 

-.5608631 49492627 D-5»_ 

4 

.1327389*0843400-6, 

4 

.9169758450660-8, 

4 

-.868388803710-9, 

4 

-.3780730050-11, 

4 

.6631455860-11 » 

4 

-.505843900-12, 

4 

-.27207820-13, 

4 

.9853310-14, 

4 

-.793980-15 , 

4 

-.67570-16, 
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* .26250-16, 

• -.2900-17, 

* -.150-18, 

• .100-18, 

* -.20-19/ 

X*8.D0/U 
6111*1.00 
GC2l«2.DO*X-l.DO 
R f 1 > » A < 1 1 ♦A(2)*G (2 ) 

SCI I*BC1)+BC2I*GC2 » 

R( 2 l*CCi)*CC2)«GC2 > 

SC2)«OCl)*DC2)«GC2> 

00 10 K*3,19 

G(3)*(4.00«X-2.00)«G(2)-G(1 ) 

R 1 1 1 «R 1 1)4A(K)*G(3) 

SC 1 l*S Cl l»BCK |*GC 3 I 
RC2)«RC2t*CCK)«GC3> 

SC2)*SC2H>DCKt»GC3) 

GC 1 ) *GC 2 ) 

G C 2 ) *G C 3 ) 

1C CONTINUE 

T*. 7978 8456C 802 8654D0/0SQRT C U> 

SN*DS I NCU-. 7853981 633974483DCI 
C S*0C0SCU-. 7853 981 63397448 3DC I 
GO T0C2C, 301, KIND 
2C BS JYC 1 )*T«CR ( 1 |«CS-S C 1 )»SN) 

BSJYC 2 l*T«CRC 2 )*SN+S C 2)*CS » 

GO TO 40 

30 BSJYCl )*T«(SC1)*CS+RC1)«SN) 

BSJYC 2 )*T*{ S C2 )*SN-R (2 l*CS » 

4C RETURN 

END 

SUBROUTINE JNSCJJ,U,M) 

INTEGER K , K A , KK « H 

DOUBLE PRECISION A , B , DC 2 > , DM ,G C 2 5 I , J J C M I • 

c P C 3 ) , 0( 3 I ,U 

DM*2*M 
PC 1 )*0.00 
0 C 1 ) *1 • DC 
PC 2 1*1 .00 
OC 2 l*DM/U 
DC1I*PC2)/0(21 
A*2 .00 

12 3*(DM+A»/U 

PC3)*B*PC2»-PCi I 
0 ( 3 ) *B*Q C 2 l-OC 1 ) 

0 C 2 )*P C 3 ) /QC 3 I 

I FC DABS CDC 1 )-DC 2 > » .LT.l .0-15 ) GO TO 20 

P C 1 l*PC2l 

PC2I-PC3I 

0C1)*QC2) 



eo 


2C 


30 

4C 

c 

c 


1C 

20 


30 


40 

50 

C 


6C 

7: 

BC 


C 

c 


0(21*0(3) 

0(1 )*0( 2 I 
A*A*2.D0 
GO TO 10 
G(M)*D(2) 
KA«H-2 

00 30 K*ltKA 
KK*M-K 


A»2*KK 

G(KIC)*U/(A-U*G(KK + 1 ) ) 

IF(G(KK).EQ«0.00) G(KK)*1.0-35 

CONTINUE 

DO 40 K»2*M 

JJ(K4ll«C(K)*JJ(K) 


CONTINUE 

RETURN 

END 

SUBROUTINE 

PURPOSES 


BSL2(S0L»U*B$IK*N> 

TO COMPUTE MOO I FI ED BESSEL 
TYPE FOR ORDERS 0 THRU N-l 


FUNCTIONS* * I * OR **' 
IN DOUBLE PRECISION. 


INTEGER N* SOL 

DOUBLE PRECISION BSIK(N!,U 
I F( S0L.E0.2 )G0 TO 30 
I F(U.GE. 8.000 ) GO TO 10 
CALL I0I1(U*BSIK) 

GO TO 20 

CALL LUKE2(U.S0L»aSIK) 
IF(N.LT.2)RETURN 
CALL INS(BSIK.UtN) 

RETURN 

IF(U.EQ.O.ODC) GO TO TO 
I FIU.GE.5.C00 ) GO TO 40 
CALL K0K1(U»BSIKI 
GO TO 50 

CALL LUKE2(U,S0L,BSIK> 
IF(N.LT. 2 (RETURN 
RECURRENCE FORMULA 
NN*N-1 


9S 1KI K*1 ) *2*ODO*OBt 6 1 FL0AT(K-lM*BSIKIK)/U*3SIK<lt*l) 

CONTINUE 

RETURN 

DO 80 K* 1 » N 

BSIK(K)*l.C0*75 

CONTINUE 

RETURN 

END 

lo^SsirTo'EUw^E’-io' »no *ii* eesssi functions 

BY SUMMING THE SERIES. 

DOUBLE PRECISION BI(2)»T(£)*X* DMA XI * DABS 


TCl)«X/2.000 
Sill )«1. 000 
B! (2 t*T( 1 1 
TC2I«TC1I**2 
TC3 1*1 *000 
TC4I*1 .000 

TC4)*TC4)«TC2)/TC3)*«2 
Bid )*B I C 1 I +T C 4 > 

TC 5)*TC4)«TC1 | / IT C 31 ♦1.0001 
BIC2)*BIC2)*TC5) 

IFf DMAX1 CDABSCTC4) I » DABS € T C 5 1)1 .LT. 1 .00-1 5 I RETURN 
TC3I«T13)+1 *000 
GO TO 10 
END 

SUBROUTINE K0K1 C X . BK ) 

PURPOSE: TO EVALUATE 'KG' AND • K1 • BESSEL FUNCTIONS 

BY SUMMING THE SERIES. 

DOUBLE PRECISION T ( 1 0 ) . X . BK ( 2 ) »0MAX1 , DABS ♦ DLOG 

TCI )«X/2.000 

T( 2 l*T( 1 1**2 

BKC 11*1.000 

BKC2 )*TC 1 I 

T C 7 I *0 .COO 

T f 10 > «-T ( 1 I 

T ( 3 ) *0.000 

T I 4 I *0 . 0 DO 

T C 5 ) * 1 .000 

TC3)*TC3H>1.0DC 

TC4)*TC4 Ul.CDO/TI 31 

TC5)*TC5)*TC2)/TC3>992 

BKC1 I = BK C 1 I+TC5) 

TC6)*TC5)*TC4> 

TC7)*TC7)+TC6) 

T(3»*TI5I»TCll/ITC3l+l.C00l 
BK C 2 I *BK C 2 ) +T C 3 ) 

TC9)*-TC81*C2.0009T(4l+l.CDC/CTC3)4i.DDOII 
T C 1 0 ) *T ( 1 0 )+TC 9 ) 

IFCDMAX1CDABSCTC6) )tDA3SCTC9)n.GE.l. CD-15) GO TO 10 
TC21*.5772156649015329D0*DLDGCTC 1 I I 
8KC1 I — 6KC1 l*T(2UT(7) 

BKC 2 I* BKC 2 )*TC 2 )*1 .ODO/X+TC 1C 1/2.0D0 

RETURN 

END 

SUBROUTINE LUKE 2 C U » K I NO , BS I K ) 

PURPOSE: TO EVALUATE MODIFIED 8ESSEL FUNCTIONS. 10 AND II 

OR KO AND K1 FROM SHIFTED CHEBYSHEV SERIES. 

DOUBLE PRECISION A C 3 4 ) , 9 ( 2 1 ) » C C 34 ) , DC 21 ) ♦ U. 3S I K C 2 ) . X , G ( 34 ) , 
C RC2)»SC2).DEXP.DS0RT 

DATA A/1. 0032792C 5458 74D0. .84451226249209430-2. 

*. 17270063077756 65 D-3* .7247591099958960-5. .5 1358772 63 7BG2D-6 
9.56816965309120-7. .851 30 91 222 8 50- 8 ..12 3842 53640-8. 
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10 


15 


*.293016723D-10»-.7895669832D-10»-.3312712763D-10, 

•-.4497338640-11,. 179979030-1 l,. 965748320-12.. 38604240-13, 

*-. 104039340-1 2 »-• 23950450-13, .9554470-14, .4443 150-14, 

*-.858640-15, -.708780-15,. 86760-16,. 111940-15, -.12110-16, 
*-.181 30-16,. 2490-17,. 2990-17, -.620-18, -.490-1 8,. 160-1 8, 

*.70-19, -.40-19,-. 10-19, .10-19/ . 

DATA B/. 9884031 742 3C 825 800, -.11310504464692 820-1, 

*.26953261276272370-3, -.11 106685196665350-4, 

*. 632575108500490-6, -.4504733764110-7,. 3792996455680 8, 
*-.364547 1792 10- 9,. 39043755760-10. -.4579936220-11, 

*.5 3 0810630-12, -.78832360-1 3,. 11360420-1 3. 

•-.1726970-14,. 275450-15, -.45 890-16.. 7960-17, -.1430 17. 

*.270-18, -.50-19, .10-19/ , 

OATA C/. 97 58006023 26285900,- .2446744296327638 0-1, 

*-.277 2 053607 6332890- 3, -.9 73 2 146728020130-5, 

•-.629724238639810-6, -.65961142154240-7, -.961387291940-8, 

*-. 140 1 1409C1030-3, -. 475631 6654D-1C.. 81530681070-10, 

*.35408148320-10 ,.5102564070- 11, -.1804409340-11, 

*-.1023594470-11, -.52677840-13.. 107094190-12,. 261 19760-1 3, 
•-.9561290-14,-. 4713350-14,. 829240-15,.742620-15,-. 80450-16, 
*-.116570- 15,. 11070-16,. 18840-16, -.2330-17,-. 3110-17, 

510-15. 

DATA 0/1 .0359 50 858772358D0, .3546529 12433 31 110-1, 
*-.46847502816688860-3,. 161 85C63810053430-4* n-8 • 

*-. 845172048123680-6,. 5713221 8102840-7,- .4645554606610 8 , 

43 5417338570- 9, -.4575729704 0-10, .5288132810-11, 

*-•662612930- 12, .8904792 0-13, -.1272 60 7 0-13, .192086 0-14, 

*-.30 45 0-15, .50450-16, -.871 0-17, .1 560-17 ,-.290-1 8 , 

*•60-19,-. 10-19/ 

IFJKJN0.EQ.21 GO TO 20 

X*8.000/U 

G(1 1*1.000 

GJ 21*2. 000*X-1. 000 

N*34 

00 10 K*3 , N 

GCK1*J4.0DC*X-2.0D0 1*G(K-1 1-GCK-2I 

CONTINUE 
RJ 1 1*0.000 
R ( 2 1*0.000 
DO 15 K* l , N 
I *N+1-K 

R ( ll*Rf 1 1*AC 1 1*GC 1 1 
R(21*R(21*C< Il*G( I 1 

9S 1KJ 1 1 1*. 3 93942 2 30 4C143270G*RC 1 t*0£XP ( U 1/DS QRT C U1 

3SIKC21*.39894223C4014327D0*RJ21*DcXP(U1/0SQRT(U1 

RETURN 

X*5.000/U 

G( 11*1.000 

GC 2 1*2 .000*X-1 .COC 

N*2 1 


20 


o o o 


63 


00 25 K«3»N 

GCKI»C4.0D0*X-2.0D0)*GCX-1 )-GCX-2) 

25 CONTINUE 

SC1I-0. ODO 
sm-o.ooo 

00 30 K«1»N 
I«NM-K 

S( 1 )«S(1 MBC I l*G( I ) 

S(2l*S(2UDin«GU) 

30 CONTINUE 

BSIXCl)*l.253314l373155D0*S( 1 )*DEXP (-U ) /DSQRT l U I 
BS IX( 2 I *i • 25331 M3 731 5500*$ ( 2 )*DEXP(-U)/D5QRT ( U) 
RETURN 
END 

SUBROUTINE INSCIltU'M) 

PURPOSE! TO EVALUATE BESSEL FUNCTIONS OF HIGHER 

OROERS BY A CONTINUED FRACTION EXPANSION 
METHODS. 

INTEGER X»XA|XX»H 

DOUBLE PRECISION A , 9, Of 2 ) t DM ,GC25 I • 1 1 1 M I ♦ P C 3 ) »Q( 3 I , U 
I F (U.EQ.O.CDO ) GO TO 50 
DM«2*M 
PCU«0.CD3 
0 ( 1 )«l.GDO 
PC2)*1.C00 

0 C 2 ) *OM/U 
DC 1 )*P( 2 )/QC 2 ) 

A-2.000 

10 B*COM*A)/U 

PC3)«3*PC2MPCl ) 

Q( 3 )«8*Q( 2 I ♦Q( 1 ) 

DC 2)*P(3 )/Q(3) 

IFCDA3SC0m-DC2H.LT.l.0D-15) GO TO 20 
PC1)«PC2) 

PC2)«PC3) 

0 f 1 > »0 1 2 » 

Q( 2 )*Q ( 3 ) 

DC 1 1*0(2 ) 

A*A*2.CD0 
GO TO 10 
20 G ( M ) *D ( 2 ) 

XA«M-2 

DO 30 X* 1 t X A 

XX*M-X 

A«2*XX 

G(XX)*U/(A+U*GCKK*1 ) ) 

IFCDABSCGCKK) I.LE.1.0D-35 ) G ( XX )*1 .00-35 
30 CONTINUE 

DO 40 X*2 » M 

IFCDABSCIICXM.LT. 1. CD-35) GO TO 35 

1 I ( K ♦ 1 )*G(X)«I I (X ) 
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GO TO 40 

35 1 1 (K+l >*0«0D0 

40 CONTINUE 

RETURN 

50 00 60 1*3 » M 

60 11111*0.000 

RETURN 
ENO 

DOUBLE PRECISION FUNCTION CERAD(QO»R»DERIV»PS* ARtNI 
C PURPOSE* TO COMPUTE A MODIFIED HATHIEU FUNCTION 

r ( OR OERIVATIVEI OF FIRST KIND CORRESPONDING 

C TO EVEN HATHIEU FUNCTION 4 CE FUNCTIONS! 

EXTERNAL C2NP»C2NlPtC2NNfC2NIN»DC2NP»DC2NIP»DC2NNtDC2NIN 

OOUBLE PRECISION ABC 25 ! »QD*FS» P SP t OUTPUT » AR C 2 5 ! 

INTEGER RiCASE.DER1V»FNC 

COMMON/ NT ERH/N1 

COMMON/ LOCAL /DUMMY ( 8 ! * AB 

N1«N 

PSP»PS 

DO 5 I *1 1 N 

ABU l*AR< I ) 

5 CONTINUE 

CASE«M00<Rt2l«-l 
IFCQD.LT.O.ODO! CASE*CASE«-2 
IFCDERIV.EQ.l! GO TO 50 
GO TOUC*20»30»40)»CASE 
C THE VALUE OF CE2N(Z»Q! 

1C CALL SIGMA<C2NPtOUTPUT) 

CERAD*PS* OUTPUT /ABfl 1 
RETURN 

C THE VALUE OF CE2N+l(ZtQ! 

20 CALL SIGMACC2NIP.0UTPUT) 

C ERA D* PS* OUTPUT /AB ( 1 ! 

RETURN 

C THE VALUE OF CE2NC Z»-Q1 

30 CALL SIGMACC2NN .OUTPUT! 

CERAD«PSP*OUTPUT/ABC 1 ) 

RETURN 

C THE VALUE OF C2SN+1CZ.-Q) 

40 CALL SIGMACC2NIN, OUTPUT! 

CERAO«PSP*OUTPUT/ABI 1 ! 

RETURN 

FOLLOWING ARE DERIVATIVES OF FUNCTIONS 
0 GO TQC60.70f80.601 »CASE 

THE VALUE OF CE2NMZ.Q! 

0 CALL S I GMAC DC2NP. OUTPUT ) 

CERAD»PS*OUTPUT/AB( 1 ) 

RETURN 

C THE VALUE OF CE2N + 1MZ.0! 

70 CALL $IGMAC0C2NIP»QUTPUT) 

CFR AD*PS*OUTPUT /AB C 1 ) 
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RETURN 

C THE VALUE OF CE2N*(Zt>0) 

80 CALL SlGMACDC2NNtOUTPUTI 

CERA0*PSP*0UTPUT/ABC1I 
RETURN 

C THE VALUE OF CE2N+1 • (Z»-Q » 

9C CALL SIGMAIDC2NIN, OUTPUT) 

CERAD»PSP*OUTPUT/AB( 1 ) 

RETURN 

ENO 

DOUBLE PRECISION FUNCTION C2NPCK) 

C PURPOSE: TO CALCULATE K TH TERM IN SUM OF SERIES 

C FOR CE2NCZ#Q). 

00U3LE PRECISION AS ( 25 1 ♦ BS JV1 » BS JV2 1 BSYV2 .03S JVl . OBS J V2 , 

* DBSVV2 
COMMON/ LOCAL/OUMMYC 8),AB 

C0MM0N/RA0IAL/BSJV1C251 # BS J V2 < 2 5 > » BS YV2 C 2 5 > »DBSJVH25I» 

C D3SJV2C25),DBSYV2C25) 

C2NP«ABm*BSJVim*BSJV2(KI 

IF<MODCK»2I.EQ.O>C2NP*-C2NP 

RETURN 

ENO 

DOUBLE PRECISION FUNCTION C2NIPCICI 
C PURPOSE: TO CALCULATE K TH TERM IN SUM OF SERIES 

C FOR C E2N + 1 ( Z »Q ) • 

DOUBLE PRECISION A3 < 25 I # BS JVl , BS JV2 t 3SYV2 .DBS JVl »DSS JV2 » 

* DBS YV2 
COMMON/ LOCAL/ DUMMY C 8 ) ? AB 

COMMON/ RAD I AL/BSJV1I 2 5 ) » BSJV2C 251 »3SYV2 1251 ,08SJV 1(25) , 

C D3SJV2I25>»DBSYV2(25I 

C2 NIP* ABC Kl*( BS JVl <K)*BSJV2 CK*1 l+BSJVKK+1 )*3SJV2(K) ) 

IF(M0DCK,2 I.EQ.9IC2NIP — C2NIP 

RETURN 

END 

DOUBLE PRECISION FUNCTION C2NNIK) 

C PURPOSE: TO CALCULATE K TH TERM IN SUM OF SERIES 

C FOR CE2N ( Z »-Q ) • 

DOUBLE PRECISION AB 1 25 > » BS I VI , BS I V2 # BSKV2 »DSS I VI f OBS I V2 t 

* DBSKV2 
COMMON/ LOCAL/ D'JMMYC 8) » AB 

COMMON/PAO IAL/BSI VI (25 ) tBSI V2C25I »3SKV2(25)t9BSIVlC25)t 

* OBSIV2C25) »D9SKV2(25) 

C2NN*AB(K»*BSIVim*3S!V2(lO 
IFCMODCKf 2 I.E0.9 )C2NN*-C2NN 
RETURN 

END 

DOUBLE PRECISION FUNCTION C2NINCK ) 

C PURPOSE: TO CALCULATE K TH TERM IN SUM OF SERIES FOR CE2N*1(Z»-Q> 

00U3LE PRECISION AB ( 25 ) ♦ BS I VI » BS I V2 » BSK V2 »D3S I VI » DBS I V2 « DBS KV2 
C OMMON/ LOCAL/ DUMMY ( 8 )» AB 

COMMON/RAO IAL/BSIVl(25)»BSIV2(25)»3SKV2(25)*DSSIVl(25)tD8SlV2(2?)» 
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C 

C 


C 

c 


c DBSKV2 ( 25 I 

C2NIN«AB(K)*(BS IV1(M*BSIV2(M1I^BSIV1(K*1I*6S1V2(IIII 

IF(N00(K»2 > .EQ*0)C2NIN«-C2NIN 

RETURN 

END 

DOUBLE PRECISION FUNCTION DC2NP(ICI 

PURPOSE* TO CALCULATE K TH TERN IN SUN OF SERIES 

OOUBLE PRECISION ABl 2 5 ) ♦ BS JV1 t BS JV2 » BSYV2 »0BS JV1 » DBSJV2 * 

* D9SYV2 »V1 1 V2 

C0MM0N/L0CAL/DUHMY1(4),V1,V2,AB Mduioit 

CONNON/RAOI AL/BSJV1 (25)»BSJV2( 25I»BSYV2I25I»DBSJV 1(251* 

* 0BSJV2(25)*DBSYV2(25I 

0C2NP*A9(M*C-O8SJVl(KI*BSJV2(KI*Vl*9SJVl(K»*OBSJV2(KI*V2l 

IF (NOOCR *2 I .E0*D)0C2NP«-0C2NP 

RETURN 

ENO 

OOUBLE PRECISION FUNCTION DC2NIP(M 
PURPOSE* TO CALCULATE K TH- TERN IN SUN OF SERIES 
FOR CEZN*1 MZ *Q >• 

DOUBLE PRECISION AB ( 25 I . BS JV 1 * BS JV2 * BS YV2 . DBS JV1 .DBS J V2 , 

« DBSYV2»V1*V2 

CQMMON/LOCAL/DUMMY1 ( 4) » VI *V2 *AB _ Deiu inn 

COMMON/RADI AL/BSJV1 (25) »BSJV2( 25) »BSYV2( 251 tDBSJVl ( 25 ) » 

* OSSJV2(25),DBSYV2(25) 

DC2NIP*AB(KI*(“D3SJV1(K)*BSJV2(K+1 )*V1* 

« BSJVinO*09SJV2(K + l)*V2- 

9 DBS JV1 (K+l |*BSJV2(K)*V1+3SJV1(K+1)*0BSJV2(K)PV2 ) 

IF(MOD(K*2).EQ*C)DC2NIP*-DC2NIP 

RETURN 

end 

OOUBLE PRECISION FUNCTION DC2NN(K) cco ,_. 

PURPOSE* TO CALCULATE K TH TERN IN SUM OF SERIES 

DOUBLE PRECISION AB ( 25 1 » BS I V 1 » BS I V2 * BSKV2 *03S I V 1 * DBS I V2 t 
« DBSKV2 » VI » V2 

CONNON/ LOCAL/OUHNY1 ( 4 ) * VI » V2 » AB 

COHMON/R AO I AL/3SIV1(25)»8SIV2(25)*3SKV2(25)»DBSIV1(25)» 

9 DBS tV2< 25 ) »D9SKV2(25) 

OC2NN» AB ( K)*(-OBSIVl(R)*BSIV2(R)*Vl+3SIVI(Kl*DBSIV2(K)*V2) 
IFINODIKf 2) .EQ*0)OC2NN«-DC2NN 

RETURN 
E NO 

nniiBLF PRECISION FUNCTION DC2NINOC) 

PURPOSE* TO CALCULATE K TH TERH IN SUM OP SERIES 
FDR C E2N+1 ' ( Z t-Q ) • 

DOUBLE PRECISION A8(25I*3SIVI*BSIV2*3SRV2*DSSIV1*DBSIV2* 

* D3SKV2,V1»V2 

* OBSIV2(25) » D3 SKV2 (251 
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DC2NINaA8(K)*(>DBSIVl(K)*BSIV2(K*l)*Vl+ 

* BSIVICK)*DBSIV2CK+1 >*V2- 

* OBS ! VI CK+l )*3SIV2( K)*V1*B$I VI CK+1 )*OBS I V2( K)*V2 » 
IF(NODCKt2).EQ.O)DC2NIN«-DC2NIN 
RETURN 
END 

OOUBLE PRECISION FUNCTION SERADC QD »R »DER I V» PS • AR *N ) 
PURPOSE* TO COMPUTE A MODIFIED MATHIEU FUNCTION 

(OR DERIVATIVE) OF FIRST KIND CORRESPONDING TO 
ODO MATHIEU FUNCTION ( SE FUNCTIONS) 

EXTERNAL S2N2P»S2NlP»S2N2NfS2NlNtDS2N2PfDS2NlPfOS2N2N» 

* OS2N1N 

DOUBLE PRECISION AB C 2 5 ) » QD» PS t PSP » OUTPUT* AR ( 2 5 ) 

INTEGER RtCASEvDERIV 
COMMON/NTERM/N1 
C OMMON/ LOCAL/ DUMMY ( 8 ) tAB 
N 1 «N 
PSPaPS 
DO 5 I »1 t N 
A B ( I I « AR ( I ) 

5 CONTINUE 

CASE«MOD(R»2)M 
IFCQD.LT. O.OOO) CASE«CASc+2 
IF(DERIV.EQ.l) GO TO 50 
GO TOl 10t20y30»4Q) .CASE 
C THE VALUE OF SE2N+2CZ»Q) 

!•" CALL S IGMA(S2N2P»OUTPUT ) 

SERAD«-PS*OUTPUT/ABC I ) 

RETURN 

C THE VALUE OF SE2NMCZ.Q) 

20 CALL SIGMACS2N1P, OUTPUT) 

SERA D*PS*OUTPUT /ABC l ) 

RETURN 

C THE VALUE OF SE2N+2C Z»-Q) 

30 CALL SIGMACS2N2Nf OUTPUT ) 

SERAD*P$P*OUTPUT/ABC I ) 

RETURN 

C THE VALUE OF S2EN«-l(Z t -Q) 

4C CALL SIGMA (S2N1N»0UT PUT) 

SERAD*PSP*OUTPUT/ABC 1 ) 

RETURN 

C FOLLOWING ARE DERIVATIVES OF FUNCTIONS 
5C GO TOC 60*70»8C » 9Q ) .CASE 

C THE VALUE OF SE2N+2* CZ»0) 

6'> CALL SIGMACDS2N2P»OUTPUT) 

SERAO*-PS* OUTPUT/ABC 1 ) 

RETURN 

C THE VALUE OF SE2N*1MZ«0) 

7C CALL SIGMACDS2NlP f OUTPUT) 

SERAO«PS*OUTPUT/ABC 1 ) 

RETURN 
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C 

80 


C 

90 


C 

c 


c 

c 


c 

c 


c 

c 


THE VALUE OF SE2N*2»U»-Q» 
CALL SIGMAIOS2N2N. OUTPUT! 
$ERAO»PSP*OUTPUT/ABI i! 
RETURN 

THE VALUE OF SE2N+1 * C2*-0» 
CALL SIGMAIDS2N1N»0UTPUT! 
SERAD*PSP* OUTPUT /ABC1! 


RETURN 
E N 0 

sun op snifs 

, 0BSYV2 

i i » *» • ;? s VV2 • 25 • * oes -*'' 1 ‘ 25 * * 

* S2N2P .«,K,,.?^^5iUSi;iU J v l .R^,. SSJ v J ,Kn 

IFIH0DCK.21.EQ.01S2N2P— S2N2P 

RETURN 

END 

^^ E r EC T s s j^c^T^:rRri N ,«■ 0E «»,« 

OOUBLE PRECISION 2 ABi25!?BSJVl,BSJV2,aSVV2.DBSJVl.OBSJV2, 

* 0BSYV2 

jassjiSiCTii! si*.:s» .«. «..«*««». . 

IF<HOD(K»2».EQ.O) S2N1P«-S2N1P 
RETURN 

OOUBLE PRECISION FUNCTION S 2 N2N(K 1 SERI C S 

purpose: to calculate k th term in sum OF SERIwS 

00U6LE P,E F c?5 I 0^ N *Ur 5 ?Us I V,..S I V 2 . E SKV 2 .0 6 S,Vl.D8S,V2. 

* DBSKV2 

sssssJtsnt ««:is «•?«;»»..««»«» .»•> . . 

% 2H2 K,.u!K!:r.ui?«!i«.««« 2 '*.*«v 2 ««' 

I F ( MOO < R» 2 1 • EQ» 0) S2N2N***S2N2N 


RETURN 

c Rjn 

double precision FUNCTION S2N1N(K| 
purpose: to calculate k th term in 

FOR SE2N+l(Z»-0). 

OOUBLE PRECISION AB I 2 5 ) » BS I VI . BS I V. 
* DBSRV2 

C OMMON/LOCAL/DU*HY1 1 8 1 » AB 


SUM OF SERIES 
,BSKV2tOBS IVlf DBSIV2» 


o n 


89 


C 

c 


c 

c 


c 

c 


CONNON/RADIAL/BSIV1I 25) »BSI V2I25) »3SKV2I25) »DBSIV1C25) » 

* D3SIV2(25)tDBSKV2(25) 
S2N1N*AB(K)*(BSIV1(K)«BSIV2(K*1 )-BSIVl(K*l)*BSIV2(K)) 

1 FI NOOIK 1 2 )«EQ*0 )S2NIN*-S2N1N 

RETURN 

END 

DOUBLE PRECISION FUNCTION OS2N2PCK) 

PURPOSE* TO CALCULATE K TH TERN IN SUN OF SERIES 
FOR S E2N + 2 * (2« Q ) • 

DOUBLE PRECISION AB( 25 ) » BS JV l * BS JV2 « BSYV2 »DBS JV 1 tDBS JV2 

* DBS YV2 » VI » V2 
CONNON/LOCAL/DUHHY1 C4)tVl »V2*AB 

CONNON/RADI AL/BS JV1 ( 25) »BSJV2(25)»BSYV2(25)»DBSJV1(25)» 

* OBSJV2C25)#DBSYV2C25l 
DS2N2P*AB(K)«(-DBSJVl(K)eBSJV2(K*2 )*V1+ 

* BSJVl(K)*D6SJV2(K+2 )*V2* 


* 03SJV1IK4-2)*5SJV2(K)*V1-8SJV1(K*2 )*D3S JV2 1 K ) *V2 ) 
I FI NOD (K, 2 ).EO.O)DS2N2P«-DS2N2P 

RETURN 

END 

DOUBLE PRECISION FUNCTION DS2N1PCK) 

PURPOSE: TO CALCULATE K TH TERN IN SUN OF SERIES 

FOR S E2N*1 ' ( 2 t Q ) • 

DOUBLE PRECISION AB C 2 5 ) » BS JV1 * BS JV2 * BSYV2 » DBS JV1 » DBS JV> 

* DBS YV2 » VI t V2 
C0HM0N/LCCAL/DUHMY1I4)»V1»V2*AB 

COMMON/RADI AL/BS J VII 25) t 8SJV2I 25) » BSYV2 C 2 5 ) t DBS JV1 1 25 ) » 

* DBSJV2I25 I »DBSYV2C25) 
DS2NIP*AB(IC)*|-DBSJV1 CK)«3S JV2 IK + I )*V1* 

* 3SJV1IK)*D3SJV2IK+1 )*V2* 

* OBSJVKK + l )*BSJV2(K)*V1-BSJV1CKM )*DBSJV2(KI*V2 I 
IFIM0DIK»2).E0.0)DS2N1P*-0S2N1P 

RETURN 


* 


» 


END 

DOUBLE PRECISION FUNCTION DS2N2NIK) 

PURPOSE: TO CALCULATE K TH TERN IN SUN OF SERIES 

FOR SE2N+2* CZf-0). 

DOUBLE PRECISION AB I 25 ) , BS 1 VI » BS I V2 t 3SKV2 .DBS I V 1 1 DBS I V2 , 

* D8SKV2 * VI t V 2 
COMNON/LOCAL/DUNMY1 U),V1 ,V2 ,AB 

COMMON/ RADI AL/BSlVl(25)»BSIV2(25lfBSlCV2f25)tDBSIVl(25)t 

* D9SIV2(25)tDBSKV2(25) 

DS2N2N«ABIK)*(-OBS IV1(K)*BSI V2IK+2 )*¥!♦ 

* BSI VI IK)*D3SlV2f K+2)*V2* 

* D3SIV1 (IU2>*3SIV2<K)*V1-9S1V1CIC+2)*03S1V?{KI*V2) 
IFIM00IK»2)*£0«0)DS2N2N*-DS2N2N 

RETURN 

END 

DOUBLE PRECISION FUNCTION DS 2N1NI K ) 

PURPOSE* TO CALCULATE K TH TERN IN SUN OF SERIES 
FOR SE2N*1M2»-Q). 



r> vh o 


90 


C 

c 

c 


5 


C 

n 


c 

20 


C 

30 


C 

40 


r 

v 


r 


DOUBLE PRECISION ABC 25 1 1 8$ I VI » BSI V2 *BSKV2*DBS I VI » DBS I V2 * 

« DBSKV2 1 VI * V2 

C0HM0N/L0CAL/DUMMYlC4!*Vl »V2*A8 

COMHON/RAOI AL/BS I VI C 25!*BSIV2C25)*BSKV2C25)»DBSIV1C25!* 

* 0BSIV2C25I»0BSKV2C25J 

DS2NIN-ABCM*C-DBSIVKKI*BSIV2CK4n*Vl* 

* 3$IV1CK1*DBSIV2CK*1>*V24 

* 0BSIVlCK4ll*BSIV2CK»*Vl-BSIVlCK4ll*0BSIV2OCI*V2» 

I F ( HOD C K » 2 I .EQ.O ! DSZNIN*-DS2N1N 

RETURN 

OOUBLE PRECISION FUNCTION FERADC QD,R .OER I V PS t AR *N I 
PURPOSE l TO COMPUTE A MODIFIED NATHIEU FUNCTION 
PURPOSE * ;U x '-“„ lvjmv5 , Qf SECONO KINO CORRESPONDING TO 

EVEN MATHIEU FUNCTION C CE FUNCTIONS! 

EXTERNAL FY2N*FV2Nl»FK2N»FK2Nl»DFV2NiOFV2Nl»OFK2N*OFK2Nl 

OOUBLE PRECISION ABI 2 51 1 0D» PS »PSP tOUTPUT* AR 1 2 5 1 1 P I 
INTEGER R t CASE* DERI V 
COMMON/ NT ERH/N1 

COMHON/LOCAl/OUMMY ( 8 ! * A 8 

DATA PI/3.141592&53589793D0/ 


Nl*N 

PSP«PS 

00 5 I*1»N 

ASCI ! * AR C I 1 

CONTINUE 

CASE*M0DCR*2 1*1 

IFCOD.LT. 0.000! CAS?*CASE+2 

IFCDERIV.EO.il GO TO 50 

GO TOC 10 * 20 * 30*40 1 »CASE 

THE VALUE OF F£Y2NCZ,Q! 

CALL SIGHACFY2N* OUTPUT ! 

F ER AD* PS* OUT PUT /ABC 1 1 
RETURN 

THE VALUE OF FEY2N*1CZ*QI 
CALL SIGMACFY2N1»0UTPUT! 
FERAD*PS*OUTPUT /ABC 1 ! 

RETURN 

THE VALUE OF FEK2NCZ*-Q1 
CALL SIGMACFK2N*0UTPUT1 
FERA0*PSP*0UTPUT/CABC11*PI 1 
RETURN 

THE VALUE OF F EK2 N+ 1 C Z »-Q ! 

CALL SIGMACFK2N1*0UTPUT1 
FERA0*PSP*0UTPUT/CABC1I*PI ! 

RETURN 

FOLLOWING ARE DERIVATIVES OF FUNCTIONS 
GO TOC 60 *70*80* 901 *CASE 
THE VALUE OF FEY2N'CZ*Q1 
CALL SIGMAC0FY2N»0UTPUT1 
FERAD*PS*OUTPUT/AB ( 1 1 


o o 


RETURN 

C THE VALUE OF FE V2N*1 * I Z » Q ) 

70 CALL S I GMA ( DFY2N1 • OUTPUT ) 

FERAD«P$*OUTPUT/AB< 1) 

RETURN 

C THE VALUE OF FEK2N' ( Zt-Q) 

80 CALL SIGMAI0FK2N, OUTPUT! 

FERAD«PSP*0UTPUT/IABI1 l*PI ! 

RETURN 

C THE VALUE OF FEK2N+1 • I Z»-Q) 

90 CALL SIGMAIDFK2N1 , OUTPUT) 

F ERAD« PSP* OUTPUT / < AB ( 1 ) *P I | 

RETURN 
ENO 

DOUBLE PRECISION FUNCTION FY2NIK ! 

PURPOSE: TO CALCULATE K TH TERN IN SUM OF SERIES 

FOR FEY2NCZ,0i. 

DOUBLE PRECISION AB ( 25 ! • BS JVI « BS JV2 • BS YV2 » DBS JVl » DBS J V? » 

* 0BSYV2 
COMMON/LOCAL/DUMMYIC 8!*AB 

COMMON/RADI AL/BS JVl ( 25 I * BS JV2C 25! * BSYV2C 25 i » 0 3S JVl ( 2 5 i * 

* 0B$JV2C25J,DSSYV2(25! 
FY2N«A3(RI«BSJV1(K!»BSYV2(K! 

IF(M0D<K,2).SQ.0)FY2N*-FY2N 
RETURN 
ENO 

DOUBLE PRECISION FUNCTION FY2N1 (K ) 

PURPOSE: TO CALCULATE K TH TERH IN SUM OF SERIES 

FEY2N + 1 ( Z t Q ) • 

DOUBLE PRECISION A B ( 2 5 ! » BS J V 1 » BS J V2 » BS YV2 » DBS J VI . OB S JV2 » 

* D8SYV2 
C OMMON/ LOCAL /DUMMY1 1 8 !» A3 

COMMON/RAOIAL/BSJVl(25!,BSJV2(251,BSYV2C25!,DBSJVl(25i, 

* OBSJV2C25! »03SYV2(25I 

FY2N1«AB(K|*(8SJV1(K1*3SYV2IK+11*BSJV1(K*1 1*BSYV2 ( K ! ! 
IFIH0D(K,2!.E0.01FY2N1«-FY2N1 
RETURN 
ENO 

D0U3LE PRECISION FUNCTION FK2NIK! 

PURPOSE: TO CALCULATE K TH TERM IN SUM OF SERIES 

FOR FEK2N(Z»Q!. 

DOUBLE PRECISION A 3 1 25 ! , BS I V 1 , BS I V2 , 9SKV2 ,OSS I V 1 , DBS I V’ , 

* 0BSKV2 
COMMON/ LOCAL/ DU MMY1I 8!» AB 

COMMON /RAO I AL/3S I VI ( 25i»BSIV2l25!*BSKV2(25!fDflSIVl(’5)t 

* 03S I V2 C 25 i »DBSKV2(25I 
FK2N«ABCK!*BStVl(Kt*BSKV2 IK) 

RETURN 
ENO 

DOUBLE PRECISION FUNCTION FK2N1IK) 

PURPOSE: TO CALCULATE K TH TERM IN SUM OF SERIES 


o o 


92 


C 

C 


c 

c 


c 

c 


DOUBLE PRECIS ION^AB tilt t BS I VI * BSI V2* B$KV2*DBS IVltOB$IV2» 

* D8SKV2 

RETURN 

END 

?ssj!)! ” , ^ $ iii!:eCSS"rt!:'?si R u * u « » f 

* DBSYV2»VlrV2 

*OMJN.»»m»l-0KJ«!”i»SYV2U»*»l*«SJVUKI*OBS¥»JUI»VJ» 
IFCMODU.2>.EQ.OIDFY2N— DFY2N 

RETURN 

END 

SSSft!.: M ?S , KSI , S N T!5'?i:i , ?i suh 0f S E «. ES 

oo U8LE M . , ers.s Y isn; , ,nsi«.« J v*..s»»*.»««‘.D.»«. 

* DBSYV2 »V1 1 V2 

« OBSJV2(25»»DBSYV2(25l 

0FY2NI.»B(*I»(-0»SJV1IKI»BSYV2IK*11«V1* 

BSJVl CK)*DBSYV2(K^1 l*V2” 

l OBSJV1Ik!i»*BSYV2CKI*V1^BSJV1IIC*1»*OBSYV2(K»*V2» 

IF1HOO(K.2».EO.O)OFY2N1 — DFY2N1 

RETURN 

END 

?^^ E : R “i S ^C^mTTr^i , ‘!N SUH OF SeRlES 

DOUBLE PRECIS10N K iBt25)?BSIVl,BS!V2,BSKV2,0BS1Vl.DBSIV2, 

* DBSKV2 »V1 f V2 

^Si^^BMJi:»;:U^» U BSKV 2 « 25 ,.OBS,Vl.: S ,. 

RETURN 
c fcJft 

;sJS!.: , ‘ i 5i s fsre!;s:5i , s N T B /s:!i , ii e «h of s E b, es 

OOUBU PR E C IS ! ON* *bI 25 W BS^Vl *8SIV2»3SKV2»0BSIVl*DflSlV2i 

. .m ^ it V 2 » V 1 # V 2 


C-2 



cr o oi o 


c 

c 

c 


5 


C 

1 : 


c 

2C 


C 

30 


C 

40 


C 0MM0N/ LOC AL/DUMMY1(4),V1,V2,AB 
CONMON/RADI AL/BSI VI 1 25 1 »BS IV2C25I f BSKV2I 251 t06SIVl(25lt 
n c K , u . 4 ,„ U( DBSIV2I25) , DBSKV2 (25 ) 

DFK2N1-ABC K)*<-D5SIVim*BSKV2CIC*n*Vl* 

* BSIV1(,c **08SKV2*K4ll*V2* 

RETURN ° BSI UK + n * BSKV2lK,0Vl " B SIVlCK4iMDBSKV2IM*V2l 
END 

PURP0SE: RE TO S rn!lpnTc CTt0N GeRA <>< QD.R , DER IV, PS ,AA ,N » 

PURPOSE: JO COMPUTE A MODIFIED MATHIEU FUNCTION 

ooo SSi?{I , ;S5 c ?IoS B i°SS ?^c°t5o°Js! spono,ng to 

oSJlLE A pRE^fsfiN V A2^5^ 2 SD , p5 2 pli D S Y2N2,OCY2N1 *° GK2N2 * DGK2N1 

C0MM0N/NTERM/N1 
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C0MM0N/L0CAL/DUMMY(8)»AB 
DATA PI/3. 1415926535897930^/ 
N1«N 
PSP*PS 


DO 5 1*1, N 
ABU I * AR { I ) 

CONTINUE 
CASE*M0D(R,2 )M 
I F( OD. LT. C • 0D0 I CASE*CASE+2 
IFIDER IV.E0.1 | GO TO 50 
GO T0( 10,20,30,40) ,CASE 
THE VALUE OF GEY2N*2(Z.Q) 
CALL S I GMA ( GY2N2 , OUTPUT I 
GERAD*-PS»OUTPUT/AB( 1 ) 
RETURN 

THE VALUE OF GEY2N*1(Z»Q) 
CALL SIGMACGY2N1, OUTPUT} 
GERAD = PS*OUTPUT/AB(l ) 

RETURN 

THE VALUE OF GEK2 N+2 ( 2, -Q } 
CALL S 1 GMA C GK2N2 , OUTPUT ) 

GER AD* P S P* OUTPUT / ( AB ( 1 ) *P I » 
RETURN 

THE VALUE OF GEK2N+1 I Z»-Q) 
CALL SIGMA(GK2N1, OUTPUT) 

GER AO* PSP* OUTPUT/ (AB(1)*PI) 
RETURN 

following are derivatives of 
GO TO( 6C ,73,80, 90 ) ,CASE 
THE VALUE OF GEY2N+2 • (Z ,Q ) 
CALL S IGMA ( DGY2N2 , OUTPUT) 
GERAD*-PS*OUTPUT/AB ( 1 ) 

RETURN 

THE VALUE OF GEY2N + 1MZ,Q) 
CALL S IGMA ( DGY2NI , OUTPUT ) 


FUNCTIONS 



o o 
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c 

80 


c 

9? 


C 

c 


c 

c 


SUM OF SERIES 


C 

C 


GER AD*PS*OUTPUT /ASCII 

?HE U VALUE OF GEK2N*2MZ,-0> 

CALL S IGHAC DGK2N2 » OUTPUT I 
gerao*psp*output/cabcii*pi> 

THE°VALUE OF GEK2N+1 * CZ»-Q> 

CALL SIGNAC DGK2N1 » OUTPUT I 

gerao«psp*output/cabci»*pi> 

RETURN 

double PRECISION FUNCTION CYEN2IK1 
PURPOSE! TO CALCULATE K TH TERN IN 

double PREC?SI0N*AB(2S^ *BSJV1*BSJV2 •BSVVZtDBS JVl *DBSJV2 t 

0 DBSYV2 

I F ( MODI R» 2 I • EQ» 01 GV2N2* 

return 

S5“E”“i”“c!iS:S"l| ,, TS Y mi‘!« SUN OP SERIES 

DOUBLE PRECISION* AaTzSItBSJVl *BSJV2»BSVV2 tOBS-IVl# DBS JV2* 

# 03SYV2 

COHHON/LOCAL/OUNMY1J8I.RB m OBSJVl(25lf 

C0HH0N/RA0IAL/BSJVH25I,BSJV,C2jI.B5YV. 

* 0BSJV2C25I.DBSYVZCZ5I , ( K , , 

%V 2 N1.»BIK..IBSJVIIM.BSY*2.R*1>-»S-<VI.NMI 

IF(H00CK»2 ),EQ.0IGY2N1«-GY2N1 
RETURN 

ga-r^KSI'^-U SUN OP SERIES 

DOUBLE PRECIS I0N K »M2M '.BS IVl.BSIV2.BSRV2>0BSIV1.0BS,V2. 

^ 0BSKV2 

ESKKSIllS&a;.;;;;'"*'”'-* 

•«„ ......msusks::..... 

RETURN 

SoUr^i'^LCUU^'rTH'TE^lN SUN OP SERIES 

DOUBLE PRECIS ,0N*AS*2M * 3S I V 1 »BS I V2 • BSKV2 • DBS I VI *06SI V2 t 

^ JU 0BSKV2 

9 


C OMMON/ LOCAL /DUMMY 11 8 )* AS 

COMNON/RADIAL/BSlVlf 25 ) » BS1 V2 ( 25 ) » BSKV2 ( 25 » « DBS I VI ( 25 ) t 

* OBSIV2(25l,DBSKV2<25» 
GK2N1*AB(K)*(BSIV1(K)*8SKV2(K+1)*BS1V1 (K + l )*8SKV2(ICM 

RETURN 

END 

DOUBLE PRECISION FUNCTION DGY2N2CIO 
PURPOSE: TO CALCULATE 1C TH TERM IN SUN OF SERIES 

FOR GEY2N+2* ( Z * Q ) • 

DOUBLE PRECISION AB C 25 > . BS J VI , BSJV2 , BSYV2 » DBS JVl » DBS JV2 » 

* DBSVV2*V1,V2 
COMMON/LOC AL/DUMM Y 1 (4)tVl*V2vAB 

COMNON/RADI AL/BSJVlf 25 I * B S JV2 ( 2 5 1 » 8SYV2 C 25 1 » DBS JVl l 25 > » 

* OBS JV2C 25) »DBSYV2( 25 I 
DGY2N2«AB< K )«{-DBS JV1 (K)*3SYV2(K*2 I*V1» 

* 3SJV1(K)*DBSYV2(K+2)*V2+ 

* DBS JV1 (K+2 )*BSYV2( K)»V1~3SJV1(K*2 )*DBSYV2(K)*V2 I 
IF(M0DCK,2).EQ.0)DGY2N2«-DGY2N2 

RETURN 

END 

DOUBLE PRECISION FUNCTION DGY2N1 ( K ) 

PURPOSE: TO CALCULATE K TH TERN IN SUN OF SERIES 

FOR GEY2N + 1 * ( 2 «Q ) • 

DOUBLE PRECISION AB( 25 I * BS JV1 » BS JV2 » 3S VV2 * DBS J VI » DBS JV2 1 

* DBSYV2 » VI » V2 

C OMMON/ LOCAL/ DUNN Y1 ( 4 ) t VI »V2 »AB 

C0HM0N/RADIAL/BSJV1I25>,BSJV2(25>,3SYV2(25),D3SJV1C25I, 

* OBSJV2C25I.D8SYV2C25I 
DGY2Nl«A3UCJ*(-DBSJVl<KI*BSYV2 (K*l l*Vl+ 

* BS JV1 m»DBSYV2(K*l |*V2 + 

* D3SJVl<K4'l|*3SYV2CK)*Vl-BSJVlfKMI*DBSYV2<IC)*V2l 
1F(H0D(K,2).EQ.0IDGY2N1»-DGY2N1 

RETURN 

END 

DOUBLE PRECISION FUNCTION DGK2N2I K > 

PURPOSE: TO CALCULATE K TH TERN IN SUN OF SERIES 

FOR GEIC2N + 2'(Z«-Q). 

DOUBLE PRECISION AB ( 2 5 I » BS I VI » 3S I V2 1 BSKV2 *DBS I VI » DBS I V2 » 

* OBSK V2 » VI # V2 
COMMON/LOCAL/ DU MHYl(4)tVltV2»AB 

COMNON/RADIAL/BSIVl(25)fBSlV2(25IfBSKV2(25IfOBSlVl(25)t 

* DBS IV2I25) »DBSKV2(25) 

DGK2N2*A8( K)*(-D3SIV1 ( K >*BSKV2 (K*2 )*V1 ♦ 

* 3SIV1(KI*DBSKV2(K*2)*V2* 

* DBS I VI (K»21*BSKV2(KJ*Vl-eSlVlCiC*2l*DSSKV2nO*V2l 
RETURN 

END 

DOUBLE PRECISION FUNCTION DGK2N1CKI 
PURPOSE: TO CALCULATE K TH TERM IN SUM OF SERIES 

FOR GEK2N*lMZ,-0>. 

DOUBLE PRECISION A B ( 2 5 I , BS I V 1 * BS I V2 1 BSK V2 .DBS I VI * DBS I V2 t 



* DBSKV2*V1»V2 
C0HM0N/LDCAl/DUMMYl<4»fVltV2tAB 

COMHON/RADI AL/BSI VH25l»BSIV2C25ltBSKV2f25)»DBSIVlt25)* 

* 08SIV2(25I*DBSAV2I25I 

0CK2N1«AB(K)*<-0BSIV1IK»*BSKV2U*1I*V1^ 

* BS IVUKI*OBSKV2CA*l l*V2- 

* DBS 1 VI (K*l l*BSKV2(KI*Vl*BSIVlllt*l l*D8SKV2C K»*V2 » 


RETURN 

ENO 
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ANISOTROPIC ELLIPTIC OPTICAL FIBERS 


Soon A. Kang, Ph. D. 
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Chicago, Illinois (1991) 

The exact characteristic equation for an anisotropic elliptic optical 
fiber is obtained for odd and even hybrid Bodes in teras of infinite 
deterainants utilizing Hathieu and modified Mathieu functions. A slaplified 
characteristic equation is obtained by applying the weakly guiding 
approxiaation such that the difference in the refractive indices of the core 
and the cladding is saall. 

The slaplified characteristic equation is used to coapute the noraalized 
guide wavelength for an elliptical fiber. When the anisotropic paraaeter is 
equal to unity, the results are coapared with the previous research and they 
are in close agreeaent. 

For a fixed value of noraalized cross-section area or aajor axis, the 
noraalized guide wavelength Wo for an anisotropic elliptic fiber is saall 
for larger the value of anisotropy. This condition indicates that aore energy 
is carried inside of the fiber. However, the geoaetry and anisotropy of the 
fiber have a smaller effect when the noraalized cross-section area is very 
saall or very large. 


